# Inference of nonlinear causal effects with GWAS summary data

Ben Dai∗ Department of Statistics, The Chinese University of Hong Kong
Chunlin Li The first two authors contribute equally to this work. School of
Statistics, The University of Minnesota  Haoran Xue  School of Statistics, The
University of Minnesota  Division of Biostatistics, The University of
Minnesota  Wei Pan  Division of Biostatistics, The University of Minnesota
Xiaotong Shen  School of Statistics, The University of Minnesota

###### Abstract

Large-scale genome-wide association studies (GWAS) have offered an exciting
opportunity to discover putative causal genes or risk factors associated with
diseases by using SNPs as instrumental variables (IVs). However, conventional
approaches assume linear causal relations partly for simplicity and partly for
the only availability of GWAS summary data. In this work, we propose a novel
model for transcriptome-wide association studies (TWAS) to incorporate
nonlinear relationships across IVs, an exposure, and an outcome, which is
robust against violations of the valid IV assumptions and permits the use of
GWAS summary data. We decouple the estimation of a marginal causal effect and
a nonlinear transformation, where the former is estimated via sliced inverse
regression and a sparse instrumental variable regression, and the latter is
estimated by a ratio-adjusted inverse regression. On this ground, we propose
an inferential procedure. An application of the proposed method to the ADNI
gene expression data and the IGAP GWAS summary data identifies 18 causal genes
associated with Alzheimer’s disease, including APOE and TOMM40, in addition to
7 other genes missed by two-stage least squares considering only linear
relationships. Our findings suggest that nonlinear modeling is required to
unleash the power of IV regression for identifying potentially nonlinear gene-
trait associations. Accompanying this paper is our Python library nl-
causal(https://github.com/nl-causal/nonlinear-causal) that implements the
proposed method.

Keywords: Transcriptome-wide and genome-wide associations, nonlinear causal
effect, sliced inverse regression, two-sample inference, invalid instruments.

##  1 Introduction

Transcriptome-wide association studies (TWAS) have successfully discovered
numerous (putative) causal genes associated with complex traits and diseases
[14], using genetic variants, typically single nucleotide polymorphisms
(SNPs), as instrumental variables (IVs) [42]. Conventional TWAS applies two-
sample two-stage least squares (2SLS) to integrate expression quantitative
trait locus (eQTL) data for gene expression and genome-wide association study
(GWAS) data for a trait of interest. For each gene, treated as an exposure,
TWAS first builds a predictive model for the expression level on the eQTL data
using the cis-SNPs around this gene as IVs. Then the predicted gene expression
is obtained with the GWAS data and then tested for association with the GWAS
trait to determine whether the gene is putatively causal to the trait. In our
motivating example, we would like to identify putative causal genes for
Alzheimer’s disease (AD) to advance its mechanistic understanding and thus its
prevention and therapeutic development.

Figure 1: Left panel. A structure plot of the proposed 2SIR nonlinear causal
model, which admits a causal effect from exposure to outcome through a
nonlinear transformation. Right panel. Estimated transformations
$\widehat{\phi}$’s of TOMM40 (a well-known AD-related gene) based on 2SLS,
PT-2SLS, and the proposed 2SIR+AIR, and the resulting $-\log_{10}(p)$-values
for the causal effect are included in the title, yielding that TOMM40 is only
identified by the proposed method. Moreover, the resulting regression $R^{2}$s
for the stage one equation on
$\widehat{\phi}(x)\sim\widehat{\bm{\theta}}^{\intercal}\bm{z}$ are 0.230
(2SLS), 0.230 (PT-2SLS), and 0.253 (2SIR+AIR), which indicates that the
nonlinear model (1) is suited for this data.

Most of the previous works in TWAS [11, 14, 47] assume a linear relationship
between cis-SNPs and gene expression in the first stage, and a linear
relationship between gene expression and a GWAS trait/outcome, ignoring
potentially nonlinear effects [24]. On the other hand, to our knowledge, none
of the existing non-parametric IV regression methods are applicable to GWAS
summary data, while individual-level GWAS data are usually unavailable due to
the privacy and logistic issues, presenting challenges to incorporating
flexible nonlinear models into TWAS with GWAS summary data. In our motivating
example, the individual-level AD GWAS data from many sub-studies are
unavailable, but its meta-analyzed summary data are available. Some recently
proposed methods [45, 31] relax the linear assumption in the first stage.
However, to our knowledge, all existing TWAS approaches assume a linear or a
parametric non-linear relationship between a gene’s expression level and a
GWAS trait/disease in the second stage, which is the focus of causal
inference. Misspecification of a nonlinear effect as a linear (or other
specific) one may distort subsequent causal inference, damping the statistical
power of the TWAS method. For illustration, we consider the eQTL data for a
well-known AD-related gene, TOMM40, from our real data example; see Section
4.1 for more details. Figure 1 provides some compelling evidence for the
nonlinear effects in both stages of TWAS. In the first stage, it displays a
nonlinear relationship between the cis-SNPs and the gene expression level of
TOMM40, as evident by the higher $R^{2}$ value of the nonlinear model over
those of its linear competitors. In the second stage, a nonlinear causal
association of TOMM40 with the AD risk is strongly corroborated by the highly
significant p-value ($<10^{-9}$) obtained with our method. Consequently, this
well-known AD gene is successfully identified by our proposed method
(2SIR+AIR) but missed by both 2SLS and its power-transformed extension
(PT-2SLS), suggesting the necessity of flexible nonlinear modeling in TWAS.

Moreover, as an IV regression method, conventional TWAS relies on three key IV
assumptions to remove the confounding effects: (IV1) the IVs are associated
with the exposure, (IV2) the IVs are not associated with the outcome
conditional on the exposure, and (IV3) the IVs are not associated with the
unmeasured confounding variables. While (IV1) is straightforward to handle,
(IV2) and (IV3) are fragile in practice due to the widespread pleiotropy of
SNPs [36]. This phenomenon refers to the situation when an SNP affects the
GWAS trait/disease not mediated through the exposure, violating (IV2) and/or
(IV3) and causing severe bias in causal inference. A line of recent works [18,
39, 13] has been focusing on the remedy to the violation of (IV2) and/or
(IV3). Of note, these methods use one-sample linear models, and their
nonlinear counterparts remain unexplored.

Besides TWAS, Mendelian Randomization (MR) is another important and popular
subject in genetics that uses SNPs as IVs to infer a causal relationship
between an exposure and an outcome, typically two complex traits [27, 41].
Both TWAS and conventional MR are two-stage IV regression methods for causal
inference and they share many similarities, yet their implementations are
different due to distinct types of data being used. Although both TWAS and
conventional MR use GWAS summary data in the second stage, in the first stage
conventional MR uses GWAS summary data of sample size typically in tens of
thousands or even larger, while TWAS typically uses individual-level data of
sample size in a few hundreds or at most one or two thousands. Usually, SNPs
being used in TWAS are around the target gene (i.e. cis-SNPs) and are
correlated, while most MR methods use independent SNPs from the whole genome.
Due to these distinctions, the existing typical MR methods do not fit the TWAS
analysis.

In a nutshell, nonlinear modeling that is robust to the violation of IV
assumptions and at the same time leverages large-scale GWAS summary data lacks
for TWAS analysis as for our motivating data. To more effectively analyze our
motivating data while addressing the limitations of existing methods, we
develop an approach with the following novel aspects.

  * •

We propose a flexible model to admit an arbitrary _unknown nonlinear_ causal
relationship between an exposure and an outcome. Importantly, the proposed
model is applicable to GWAS summary data while being robust to invalid IVs.

  * •

On the ground of the proposed model, we decouple the estimation of a _marginal
causal effect_ and a _nonlinear transformation_. To assess the marginal causal
effect, the estimation and inference procedures are established based on the
proposed two-stage instrumental regression (2SIR) without estimating the
transformation. Specifically, in the first stage, sliced inverse regression
(SIR) is committed to identifying the effects of SNPs on gene expression, and
in the second stage, a sparse regression is conducted for simultaneous invalid
IV selection and unbiased estimation of the marginal causal effect. To assess
the average treatment effect, the unknown nonlinear transformation can be
consistently estimated by the proposed adjusted inverse regression (AIR).

  * •

Our theoretical result guarantees the validity of hypothesis testing and
interval estimation of the proposed method. An extensive simulation study
confirms our theory.

  * •

An application to our motivating data, the ADNI (individual-level) gene
expression data, and the IGAP AD GWAS summary data confirms both the necessity
and power of our approach. The results (Section 4.1 and Figure 5) suggest that
the proposed method can replicate the significant genes of AD identified by
2SLS, including the well-known APOE gene. Moreover, seven additional
(putative) causal genes are only identified by the proposed method, including
the well-known gene TOMM40, and a comprehensive literature search provides
strong evidence supporting these new findings. Our real data analysis suggests
potential nonlinearity in gene-trait associations, and nonlinear modeling is
suited to fully unleash the power of TWAS.

##  2 Nonlinear modeling of TWAS data

We denote a vector of instrumental variables as $\bm{z}\in\mathbb{R}^{p}$, a
scalar exposure as $x\in\mathbb{R}$, and a scalar outcome as $y\in\mathbb{R}$.
For example, in our case study (cf. Section 4.1), SNPs are used as
instrumental variables for a genetic exposure to identify its causal
association with the AD risk. Without loss of generality, we assume
$(\bm{z},x,y)$ have mean zero. Suppose $(\bm{z},x,y)$ satisfy a nonlinear
transformation model

$\phi(x)=\bm{z}^{\intercal}\bm{\theta}+w,\qquad
y=\beta\phi(x)+\bm{z}^{\intercal}\bm{\alpha}+\varepsilon,$ (1)

where $(w,\varepsilon)$ are the error terms independent of the instruments
$\bm{z}$, and $\beta\in\mathbb{R}$, $\bm{\alpha}\in\mathbb{R}^{p}$,
$\bm{\theta}\in\mathbb{R}^{p}$ are unknown parameters, and $\phi(\cdot)$ is an
unknown transformation.

In (1), $\beta\phi(\cdot)$ represents the influence of different exposure
levels on the outcome, which is one of our primary focus, while $\bm{\alpha}$
and $\bm{\theta}$ are nuisance parameters. In particular,
$\bm{\alpha}\neq\bm{0}$ indicates the violation of the second and/or third IV
assumptions. Generally, the effect $\beta\phi(\cdot)$ may not be identifiable
with the presence of invalid IVs. In literature, additional structural
constraints are imposed to avoid this issue. For example, if
$\|\bm{\alpha}\|_{0}<p/2$ is known a priori, then $\beta\phi(\cdot)$ becomes
well-defined [19]. Furthermore, note that $\beta$ and $\phi$ are only
identifiable up to a multiplicative scalar, even if $\beta\phi(\cdot)$ is
well-defined in (1). Thus, we fix $\|\bm{\theta}\|_{2}=1$ and $\beta\geq 0$ in
the subsequent discussion so that $\beta$ and $\phi$ are identifiable.

The benefits of the proposed model in (1) are twofold. First, (1) admits an
arbitrary nonlinear transformation $\phi(\cdot)$ across $\bm{z}$, $x$ and $y$,
relaxing the linearity assumption in the standard TWAS/2SLS. Second, it
includes the classical two-stage least-squares (2SLS; [19]) and Yeo-Johnson
power transformation 2SLS (PT-2SLS; [43]) as special cases. Specifically, 2SLS
assumes a linear relationship $\phi(x)=x$, and PT-2SLS assumes that $\phi(x)$
is normally distributed, and there exists a power transformation
$\psi_{\lambda}(x)$, c.f., (2.1) in [43], such that
$\phi(x)=\psi_{\lambda}(x)$. It is worth to mention that the proposed method
yields competitive performance against 2SLS when the linear assumption or
normality assumption holds, see Section 3. Overall, the proposed model (1) is
a natural extension of the two-stage linear model.

On this ground, Definition 1 summarizes the quantities of interest.

######  Definition 1 (Causal effect and transformation).

In (1), let $\|\bm{\theta}\|_{2}=1$ and $\beta\geq 0$. Then,
(i) $\beta$ is called the marginal causal effect;
(ii) $\phi(\cdot)$ is called the nonlinear transformation;
(iii) $\beta\phi(\cdot)$ is called the nonlinear effect function.

Specifically, $\beta$ summarizes the marginal effect of the causal influence
of the exposure $x$ on the outcome $y$, in that $\beta>0$ indicates the
presence of the causal relation, and the corresponding hypothesis testing and
confidence interval are developed in Sections 2.1 and 2.3. If the model (1) is
well-specified, the nonlinear effect function $\beta\phi(\cdot)$ in (iii) can
be used to measure the average treatment effect (ATE) between two
exposure/treatment levels. For example, ATE for (1) is defined as:

$\operatorname{E}(y\mid do(x=x_{1}))-\operatorname{E}(y\mid
do(x=x_{2}))=\beta\big{(}\phi(x_{1})-\phi(x_{2})\big{)},\quad\text{for any
}x_{1},x_{2}\in\mathcal{X},$

where $\mathcal{X}$ is a pre-given treatment range of interest and $do(\cdot)$
is the do-operator [33]. In our case study, $\beta>0$ indicates the presence
of the causal influence of a gene on the AD risk, and if the model (1) is
well-specified, $\phi(\cdot)$ represents the potentially nonlinear pattern of
a putative causal association. Both are demonstrated in Section 4.1, where the
significant genes with $\beta>0$ are displayed in Figure 5, and the estimated
nonlinear transformations for representative genes are illustrated in Figures
1 and 7.

Let $(\bm{Z}_{\nu},\bm{X}_{\nu},\bm{Y}_{\nu})$ be $n_{\nu}\times(p+2)$ matrix,
where each row $(\bm{z}_{\nu i},x_{\nu i},y_{\nu i})$, $1\leq i\leq n_{\nu}$,
$\nu=1,2$, represents an independent observation from (1). In what follows,
assume that we have two independent samples
$\mathcal{D}_{1}=\\{\bm{Z}_{1},\bm{X}_{1}\\}$ and
$\mathcal{D}_{2}=\\{n_{2}^{-1}\bm{Z}_{2}^{\intercal}\bm{Z}_{2},n_{2}^{-1}\bm{Z}_{2}^{\intercal}\bm{Y}_{2},n_{2}^{-1}\bm{Y}_{2}^{\intercal}\bm{Y}_{2}\\}$
from (1). Without loss of generality, we assume that $\bm{Y}_{2}$ is pre-
normalized as $n_{2}^{-1}\bm{Y}^{\intercal}_{2}\bm{Y}_{2}=1$. Importantly, we
require neither that all variables $(\bm{z},x,y)$ are observed simultaneously,
nor the availability of individual-level data
$(\bm{Z}_{2},\bm{X}_{2},\bm{Y}_{2})$, allowing the application to summary
statistics, like the GWAS summary data, for the second sample.

Our goal is to infer $\beta$ and $\beta\phi(\cdot)$ from the observed data
$\mathcal{D}_{1},\mathcal{D}_{2}$. In the sequel, we propose estimating the
marginal causal effect $\beta$ and the nonlinear transformation $\phi$
separately.

###  2.1 Estimation of marginal causal effect

The proposed procedure for estimating $\beta$ consists of two stages. In the
first stage, note that
$x\mathrel{\text{\scalebox{1.07}{$\perp\mkern-10.0mu\perp$}}}\bm{z}\mid\bm{z}^{\intercal}\bm{\theta}$
in (1), which coincides with a semiparametric transformation/single index
model [7, 4], and the sliced inverse regression (SIR) [21] can be used to
estimate $\bm{\theta}$ as the so-called sufficient dimension reduction (SDR)
direction. Specifically, given the dataset $\mathcal{D}_{1}$, SIR divides the
range of $x_{i}$ into $S$ non-overlapping slices
$\text{Slice}_{s}(s=1,\cdots,S)$, SIR estimates ${\bm{\theta}}$ as the
eigenvector of $\widehat{\bm{\Sigma}}^{-1}\widehat{\bm{\Gamma}}$ associated
with the largest eigenvalue:

$\widehat{\bm{\theta}}=\operatorname*{\arg\max}_{\bm{\theta}\in\mathbb{R}^{p}}\bm{\theta}^{\intercal}\widehat{\bm{\Gamma}}\bm{\theta},\quad\text{s.t.}\quad\bm{\theta}^{\intercal}\widehat{\bm{\Sigma}}\bm{\theta}=1,$
(2)

where $\widehat{\bm{\Sigma}}$ is the sample covariance matrix of $\bm{z}$, and
$\widehat{\bm{\Gamma}}$ is the between slice covariance matrix defined by

$\widehat{\bm{\Sigma}}=\frac{1}{n_{1}}\sum_{i=1}^{n_{1}}\bm{z}_{1i}\bm{z}_{1i}^{\intercal},\quad\widehat{\bm{\Gamma}}=\sum_{s=1}^{S}\frac{n_{1s}}{n_{1}}\bar{\bm{z}}_{(s)}\bar{\bm{z}}^{\intercal}_{(s)},\quad\bar{\bm{z}}_{(s)}=\frac{1}{n_{1s}}\sum_{x_{i}\in\text{Slice}_{s}}\bm{z}_{1i},$

with $n_{1s}$ being the number of samples in the $s$-th slice
$\text{Slice}_{s}$.

In the second stage, we estimate $\beta$ via a sparse instrumental variable
regression using the data $\mathcal{D}_{2}$. Specifically, note that the
second equation in (1) can be rewritten as

$y=\bm{z}^{\intercal}\bm{\theta}\beta+\bm{z}^{\intercal}\bm{\alpha}+e,\qquad
e=w\beta+\varepsilon,\qquad\operatorname{E}(e)=0,\qquad\operatorname{E}(e^{2})=\sigma_{e}^{2}.$
(3)

Recall that $\alpha_{j}\neq 0$ indicate $z_{j}$ violates the IV assumptions
(IV2) and/or (IV3). Motivated by [41], we separate the potential bias due to
invalid IVs from the marginal causal effect $\beta$ via the following sparse
regression:

$\begin{split}\min_{\bm{\alpha},\beta}\
(\widehat{\bm{\theta}}\beta+\bm{\alpha})^{\intercal}\bm{Z}_{2}^{\intercal}\bm{Z}_{2}(\widehat{\bm{\theta}}\beta+\bm{\alpha})-2\bm{Y}_{2}^{\intercal}\bm{Z}_{2}(\widehat{\bm{\theta}}\beta+\bm{\alpha})\quad\text{s.t.}\quad\|\bm{\alpha}\|_{0}\leq
K,\end{split}$ (4)

where $\|\bm{\alpha}\|_{0}=\sum_{j=1}^{p}\operatorname{I}(\alpha_{j}\neq 0)$
and $K\geq 0$ is an integer tuning parameter indicating the number of invalid
IVs. For implementation, $\|\cdot\|_{0}$ penalty can be replaced by a
sparsity-inducing surrogate penalty, such as SCAD [9], TLP [35], and MCP [44].
In our data analysis, we use the SCAD [9] as a computational surrogate; see
Appendix A.1 for details.

Taken together, the proposed procedure consists of the estimation of
$\bm{\theta}$ via a Sliced Inverse Regression, and that of $\beta$ via a
Sparse Instrumental Regression. This methodology is named Two-Stage
Instrumental Regression (2SIR), as summarized in Algorithm 1.

Input : Datasets $\mathcal{D}_{1}$ and $\mathcal{D}_{2}$

Output : Estimated causal effect $\widehat{\beta}$ and $\widehat{\bm{\theta}}$

1 (Stage 1: Sliced inverse regression) Estimate $\widehat{\bm{\theta}}$ via
(2) with $\mathcal{D}_{1}$ ;

2 (Stage 2: Sparse instrumental regression) Estimate $\widehat{\beta}$ via (4)
with $\mathcal{D}_{2}$ and $\widehat{\bm{\theta}}$;

3 (Sign adjustment for identifiability)
$\widehat{\bm{\theta}}\leftarrow\operatorname{sign}(\widehat{\beta})\widehat{\bm{\theta}}$,
$\widehat{\beta}\leftarrow|\widehat{\beta}|$;

return _$(\widehat{\beta},\widehat{\bm{\theta}})$_

Algorithm 1 Two-stage instrumental regression (2SIR) for $\beta$ estimation

Section 3 indicates that 2SIR can provide a stable and (asymptotically)
unbiased estimator in various simulated examples. Yet, in practice,
visualization of $\phi$ may shed light on the specific relationship between
the exposure and outcome. In the next section, we develop an algorithm to
estimate the nonlinear transformation $\phi$.

###  2.2 Estimation of nonlinear transformation

The challenge of estimating $\phi(\cdot)$ is twofold. First, individual-level
data of $(z,x,y)$ are usually unavailable (in both $\mathcal{D}_{1}$ and
$\mathcal{D}_{2}$), preventing the estimation of $\phi$ from the second
equation of (1). Second, $w$ is correlated with $\phi(x)$ in (1), rendering a
biased estimator when for example a least-squares regression of
$\bm{z}^{\intercal}\bm{\theta}$ is conducted over $x$.

To address these issues, we propose an Adjusted Inverse Regression (AIR) for
consistent estimation of $\phi$. An important observation is made in
Proposition 1, showing that the transformation $\phi$ is proportional to the
least-squares estimator.

######  Proposition 1.

Suppose $\operatorname{E}(\bm{z}^{\intercal}\bm{\theta}\mid
x)=\operatorname{E}(\bm{z}^{\intercal}\bm{\theta}\mid\phi(x))$ and
$(\bm{z}^{\intercal}\bm{\theta},w)$ has an elliptically symmetric
distribution. Then there exists a constant $\rho$ such that
$\phi(x)=\rho\operatorname{E}(\bm{z}^{\intercal}\bm{\theta}\mid x)$.

In light of Proposition 1, $\phi$ can be estimated by a two-stage procedure.
First, we estimate the conditional mean
$\operatorname{E}(\bm{z}^{\intercal}\bm{\theta}\mid x)$ via the least-squares
regression:

$\widehat{m}=\operatorname*{\arg\min}_{m\in\mathcal{F}}\
\frac{1}{2n_{1}}\sum_{i=1}^{n_{1}}\left(\bm{z}_{1i}^{\intercal}\widehat{\bm{\theta}}-m(x_{1i})\right)^{2},$
(5)

where $\mathcal{F}$ is a class of candidate functions. Note that (5) includes
many methods, such as smoothing spline regression [38], local regression [28],
and gradient boosting regression [10].

Then, $\widehat{\rho}$ is estimated base on the uncorrelatedness between
$\bm{z}^{\intercal}\bm{\theta}$ and $w$, that is,

$\frac{1}{n_{1}}\sum_{i=1}^{n_{1}}(\bm{z}_{1i}^{\intercal}\widehat{\bm{\theta}})\left(\bm{z}_{1i}^{\intercal}\widehat{\bm{\theta}}-\widehat{\rho}\widehat{m}(x_{1i})\right)=0,\qquad\widehat{\rho}=\frac{\widehat{\bm{\theta}}^{\intercal}\sum_{i=1}^{n_{1}}(\bm{z}_{1i}\bm{z}_{1i}^{\intercal})\widehat{\bm{\theta}}}{\widehat{\bm{\theta}}^{\intercal}\sum_{i=1}^{n_{1}}\widehat{m}(x_{1i})\bm{z}_{1i}}.$
(6)

Finally, the AIR estimator is $\widehat{\phi}=\widehat{\rho}\widehat{m}$. The
procedure is summarized in Algorithm 2.

Input : Dataset $\mathcal{D}_{1}$ and estimated $\widehat{\bm{\theta}}$

Output : Estimated $\widehat{\phi}$

1 (Stage 1: Inverse regression) Estimate condition mean $\widehat{m}(\cdot)$
via (5);

2 (Stage 2: Adjustment ratio) Compute $\widehat{\rho}$ by (6);

return _$\widehat{\phi}(\cdot)=\widehat{\rho}\widehat{m}(\cdot)$_

Algorithm 2 Adjusted inverse regression (AIR) for $\phi(\cdot)$ estimation

It is worth noting that the estimation of $\phi$, as well as $\bm{\theta}$, is
solely based on the dataset $\mathcal{D}_{1}$. Moreover, AIR allows the
estimation of a non-invertible transformation $\phi$, provided that $\phi(x)$
is a sufficient statistic for conditional mean
$\operatorname{E}(\bm{z}^{\intercal}\bm{\theta}\mid x)$. This is in contrast
to the existing literature on data transformation (see [43] and references
therein), where only invertible transformations are considered. In Section 3,
the numerical results demonstrate the advantages of our method in detecting a
quadratic relationship.

###  2.3 Hypothesis testing and confidence interval

This section presents inferential procedures for the marginal causal effect
$\beta$, including hypothesis testing and confidence intervals. Before
proceeding, Theorem 1 summarizes the asymptotic properties of the 2SIR
estimator.

######  Theorem 1.

Let $\widehat{\beta}$ be the 2SIR estimator produced by Algorithm 1 with
$\|\cdot\|_{0}$ penalty, SCAD, TLP, or MCP being used in (4). Assume
Conditions B.1 and B.2 in Appendix B. If $K=|A|$ in (4) and $(w,\varepsilon)$
is normally distributed, then

$\begin{split}&n_{2}^{1/2}(\widehat{\beta}-\beta)=|n_{2}^{1/2}\beta+\zeta-\eta|-n_{2}^{1/2}\beta+o_{p}(1),\quad\zeta\mathrel{\text{\scalebox{1.07}{$\perp\mkern-10.0mu\perp$}}}\eta,\\\
&\zeta\sim
N(0,\Omega_{X}{\sigma^{2}_{e}}),\quad\eta\sim\sqrt{r}\beta\Omega_{X}\bm{\theta}^{\intercal}\widetilde{\bm{\Sigma}}\bm{\xi},\end{split}$
(7)

where $A=\\{j:\alpha_{j}\neq 0\\}$, $n_{2}/n_{1}\to r$ and
$n_{1}^{1/2}(\widehat{\bm{\theta}}-\bm{\theta})\overset{d}{\longrightarrow}\bm{\xi}$,
$\widetilde{\bm{\Sigma}}=\bm{\Sigma}-\bm{\Sigma}_{*A}\bm{\Sigma}_{AA}^{-1}\bm{\Sigma}_{A*}$,
$\Omega_{X}=(\bm{\theta}^{\intercal}\widetilde{\bm{\Sigma}}\bm{\theta})^{-1}$,
and $\bm{\Sigma}_{*A},\bm{\Sigma}_{A*}$ denote the columns and rows of
$\bm{\Sigma}$ indexed by $A$, respectively.

Now, we infer $\beta$ based on Theorem 1. First, consider the hypotheses:

$H_{0}:\beta=0\quad\text{versus}\quad H_{a}:\beta>0,$

where rejecting the null hypothesis $H_{0}$ indicates an evidence for causal
influence of the exposure $x$ on the outcome $y$. Define the pivotal test
statistic

$\widehat{T}=\frac{n_{2}^{1/2}\widehat{\beta}}{\widehat{\sigma}_{e}(\widehat{\bm{\theta}}^{\intercal}\widehat{\bm{\Sigma}}\widehat{\bm{\theta}}-\widehat{\bm{\theta}}^{\intercal}\widehat{\bm{\Sigma}}_{*A}(\widehat{\bm{\Sigma}}_{AA})^{-1}\widehat{\bm{\Sigma}}_{A*}\widehat{\bm{\theta}})^{1/2}}.$
(8)

Given a significance level $\alpha\in(0,1)$, the null hypothesis $H_{0}$ is
rejected if and only if $\widehat{T}>\Phi^{-1}_{N(0,1)}(1-\alpha/2)$, where
$\Phi^{-1}_{N(0,1)}(\cdot)$ denotes the quantile function of $N(0,1)$.

As a consequence of Theorem 1, Corollary 1 justifies the proposed test.

######  Corollary 1.

Assume the conditions in Theorem 1. The following statements are true.
(i) Under the null hypothesis $H_{0}:\beta=0$, we have

$\limsup_{n_{2}\to\infty}P_{H_{0}}\Big{(}\widehat{T}>\Phi^{-1}_{N(0,1)}(1-\alpha/2)\Big{)}\leq\alpha.$
(9)

(ii) Under the alternative hypothesis $H_{a}:\beta=n_{2}^{-1/2}h$, we have

$\liminf_{n_{2}\to\infty}P_{H_{a}}\Big{(}\widehat{T}>\Phi^{-1}_{N(0,1)}(1-\alpha/2)\Big{)}\geq
P\left(|N(\Omega_{X}^{-1/2}\sigma_{e}^{-1}h,1)|>\Phi^{-1}_{N(0,1)}(1-\alpha/2)\right).$
(10)

Empirically, Section 3 shows that the proposed test can control the Type I
error under the null hypothesis $H_{0}$ while possessing desirable power under
$H_{a}$. Moreover, we develop a combined test over a different number of
slices for 2SIR, see Appendix A.3.

Next, we consider constructing a valid confidence interval for $\beta$.
Indeed, this can be challenging, since the asymptotics of the SIR estimator
depends on an unknown distribution $\bm{z}\mid x$ [46], which is intractable.
To overcome this difficulty, we propose a resampling-based confidence interval
in light of Theorem 1. Specifically, by the triangle inequality,
$n_{2}^{1/2}|\widehat{\beta}-\beta|\leq|\zeta-\eta|+o_{p}(1)$, and hence

$\liminf_{n_{2}\to\infty}P\left(n^{1/2}_{2}|\widehat{\beta}-\beta|\leq
Q(1-\alpha)\right)\geq 1-\alpha,$

where $Q(\cdot)$ denotes the quantile function of $|\zeta-\eta|$. For
implementation, we first compute $(\widehat{\bm{\theta}},\widehat{\beta})$ via
Algorithm 1,
$\widehat{\Omega}_{X}=(\widehat{\bm{\theta}}^{\intercal}(\widehat{\bm{\Sigma}}-\widehat{\bm{\Sigma}}_{*A}\widehat{\bm{\Sigma}}_{AA}^{-1}\widehat{\bm{\Sigma}}_{A*})\widehat{\bm{\theta}})^{-1}$,
and
$\widehat{\sigma}_{e}^{2}=n_{2}^{-1}(\bm{Y}_{2}^{\intercal}\bm{Y}_{2}-\bm{Y}_{2}^{\intercal}\bm{Z}_{2}(\bm{Z}_{2}\bm{Z}_{2})^{-1}\bm{Z}_{2}^{\intercal}\bm{Y}_{2})$.
Then the bootstrap estimates $\widehat{\bm{\theta}}^{*}_{l}$s are computed as
$\widehat{\bm{\theta}}^{*}_{l}=\operatorname{sign}(\widehat{\bm{\theta}}^{\intercal}\widetilde{\bm{\theta}}^{*}_{l})\widetilde{\bm{\theta}}^{*}_{l}$,
where $\widetilde{\bm{\theta}}^{*}_{l}$ is computed via Step 1 (SIR) in
Algorithm 1 based on resampling $\mathcal{D}_{1}$, and $\zeta^{*}_{l}\sim
N(0,\widehat{\Omega}_{X}\widehat{\sigma}_{e}^{2})$ is generated according to
its asymptotic distribution; $l=1,\ldots,M$, where $M$ is the Monte-Carlo
size. In this way, we approximate the distribution of $\eta$ by the Monte-
Carlo sample:

$\eta^{*}_{l}=\frac{1}{2}n_{2}^{1/2}\widehat{\beta}\widehat{\Omega}_{X}\Big{(}(\widehat{\bm{\theta}_{l}^{*}})^{\intercal}(\widehat{\bm{\Sigma}}-\widehat{\bm{\Sigma}}_{*A}\widehat{\bm{\Sigma}}_{AA}^{-1}\widehat{\bm{\Sigma}}_{A*})\widehat{\bm{\theta}}^{*}_{l}-\widehat{\bm{\theta}}^{\intercal}(\widehat{\bm{\Sigma}}-\widehat{\bm{\Sigma}}_{*A}\widehat{\bm{\Sigma}}_{AA}^{-1}\widehat{\bm{\Sigma}}_{A*})\widehat{\bm{\theta}}\Big{)};\quad
l=1,\ldots,M.$

As a result, the $(1-\alpha)$-confidence interval is given by

$\beta\in\left[\max(0,\widehat{\beta}-n_{2}^{-1/2}\widehat{Q}^{*}(1-\alpha)),\
\widehat{\beta}+n_{2}^{-1/2}\widehat{Q}^{*}(1-\alpha)\right],$ (11)

where $\widehat{Q}^{*}(\cdot)$ is the quantile function of
$|\zeta^{*}_{l}-\eta_{l}^{*}|$; $l=1,\ldots,M$.

It is worth noting that if $\beta=0$, then $|\zeta-\eta|=|\zeta|$ and
$n^{1/2}_{2}\widehat{\beta}=|\zeta|+o_{p}(1)$ by Theorem 1, and the test based
on (8) is thus optimal. However, in the case of $\beta\approx 0$, the
bootstrap quantile $\widehat{Q}^{*}(1-\alpha)$ is usually larger than the
quantile $Q(1-\alpha)$ of $|\zeta-\eta|$, due to additional variations of
$\eta^{*}$. Hence, the confidence interval in (11) is less efficient than the
test (9) when detecting a small signal $\beta\approx 0$.

Finally, for interval estimation of $\beta\phi(\cdot)$, there is ample
literature devoted to constructing nonparametric confidence bands in a
regression setup; see [15] and references therein.

###  2.4 Robustness to misspecified nonlinearity

The proposed nonlinear model (1) considerably relaxes the linearity assumption
in conventional TWAS methods. Nevertheless, it is possible that the nonlinear
transformation $\phi(\cdot)$ in (1) could be misspecified in practice,
especially when two structural equations do not share the same transformation
for the exposure:

$\phi(x)=\bm{z}^{\intercal}\bm{\theta}+w,\quad
y=\beta\psi(x)+\bm{z}^{\intercal}\bm{\alpha}+\varepsilon,$ (12)

where $\phi\neq\psi$ are two different nonlinear functions. In this situation,
it is generally impossible to consistently estimate $\psi$ or ATE from the
available summary statistics in TWAS. Yet, hypothesis testing in Section 2.3
remains valid.

######  Corollary 2.

In (12), if the conditions in Theorem 1 are met, then under the null
hypothesis $H_{0}:\beta=0$, we have

$\limsup_{n_{2}\to\infty}P_{H_{0}}\Big{(}\widehat{T}>\Phi^{-1}_{N(0,1)}(1-\alpha/2)\Big{)}\leq\alpha.$
(13)

As a result, in our TWAS analysis, the p-values of the putative causal genes
produced by 2SIR remain reliable regardless of whether the transformations are
correctly specified. Furthermore, the simulation results indicate that the
proposed test enables control of the Type I error and outperforms its
competitors in power in the misspecified cases; see Section 3.1 Example 6.

##  3 Simulations

This section examines numerically the performance of the proposed 2SIR and AIR
methods. Moreover, for hypothesis testing, we propose to combine tests based
on different slices, denoted as Comb-2SIR. Let $\mathcal{S}$ be a collection
of candidate slices, we combine $p$-values based on different slices
$S\in\mathcal{S}$ using the Cauchy combining method [22]:

$p_{*}=0.5-(\arctan t_{0})/\pi,\quad
t_{0}=\sum_{S\in\mathcal{S}}w_{i}\tan\Big{(}\big{(}0.5-P_{W}(|W|>\widehat{T}_{S})\big{)}\pi\Big{)},$
(14)

where the weights $w_{i}$s are nonnegative and
$\sum_{i=1}^{|\mathcal{S}|}w_{i}=1$, $\widehat{T}_{S}$ is the test statistic
in (8) with the subscript indicating its dependence on $S$, and $W\sim N(0,1)$
is a standard normal variable independent of the data. In our numerical
examples, we focus on a combined test with $w_{i}=1/|\mathcal{S}|$ and
$\mathcal{S}=\\{2,3,5,10\\}$. More discussion about the Cauchy combining
version of 2SIR over the number of slices is included in Appendix A.3.

The performance for both $\beta$ and $\phi(\cdot)$ are considered.
Specifically, the results are compared against two-stage least-squares and
two-stage least-squares based on the Yeo-Johnson power transformation (a
generalized Box-Cox transformation [43]), denoted as 2SLS and PT-2SLS,
respectively. For PT-2SLS, the optimal parameter $\lambda$ for minimizing
skewness is estimated using maximum likelihood, c.f., Section 3 in [43].

###  3.1 Testing and confidence interval for marginal causal effect

This subsection examines the proposed method for its capability of controlling
Type I and II errors, coverage, and effectiveness of confidence intervals for
simulated examples.

The simulated data $\mathcal{D}=(\bm{z}_{i},x_{i},y_{i})_{i=1,\cdots,n}$ is
generated as follows. First, $\bm{z}_{i}$ is generated independently from
$N(\bm{0}_{p},\bm{\Sigma})$, and $w_{i}=u^{2}_{i}+\gamma_{i}$, where $u_{i}$
and $\gamma_{i}$ are independently generated from $N(0,1)$. Second, $x_{i}$ is
generated as $x_{i}=\phi^{-1}(\bm{\theta}^{\intercal}\bm{z}_{i}+w_{i})$ when
$\phi$ is invertible, and $x_{i}$ is randomly selected from the solution set
$\\{x:\phi(x)=\bm{\theta}^{\intercal}\bm{z}_{i}+w_{i}\\}$ when $\phi$ is non-
invertible. Third, $y_{i}=\beta\phi(x_{i})+\varepsilon_{i}$, where
$\varepsilon_{i}=u_{i}+\zeta_{i}$, and $\zeta_{i}\sim N(0,1)$, thus $u_{i}$
acts as a confounder, and $w_{i}$ is dependent with $\varepsilon_{i}$.
Finally, the first half of individual data is provided as $\mathcal{D}_{1}$,
and the summary data $\mathcal{D}_{2}$ is produced by the second half of data
to mimic the GWAS data. Six transformations are considered: (1) linear:
$\phi(x)=x$; (2) logarithm: $\phi(x)=\log(x)$; (3) cube root:
$\phi(x)=x^{1/3}$; (4) inverse: $\phi(x)=1/x$; (5) piecewise linear:
$\phi(x)=xI(x\leq 0)+0.5xI(x>0)$; (6) quadratic: $\phi(x)=x^{2}$.

For Type I error and power analysis, we compute the proportions of rejecting
the null hypothesis out of 1,000 simulations under $H_{0}$ and out of 100
simulations under $H_{a}$, respectively. For constructing the confidence
interval, we report the averaged coverage and CI length out of 1,000
simulations. Note that the CIs for 2SLS and PT-2SLS are generated based on the
asymptotic variance in [16], the CI for the proposed method is generated based
on (11), and all CIs are left truncated at 0 since $\beta\geq 0$.

Example 1 (Standard setting). In this example, we examine the proposed method
under a standard setting. Specifically, we set $\bm{\Sigma}=\bm{I}_{p}$,
$\bm{\theta}\sim N(\bm{0},\bm{I}_{p})$ and normalize it by its norm. We
examine four cases: (i) $\beta=0$, (ii) $\beta=.05$, (iii) $\beta=.10$, (iv)
$\beta=.15$. Note that case (i) is for Type I error analysis, while $\beta>0$
in (ii) - (iv), suggesting power analysis. From (ii) to (iv), the true causal
effect is increasing, thus the power is expected to go up. Moreover, the
confidence interval is produced based on (ii) $\beta=0.05$. All empirical
results are summarized in Figure 2 (testing) and Table 1 (CI).

Example 2 (Invalid IVs). In this example, we examine the performance and
stability of the proposed method with invalid IVs. Specifically, $\bm{z}_{i}$
is generated with $\Sigma_{ij}=\nu^{|i-j|}$; $i,j=1,\cdots,p$. Then, $x_{i}$
is generated based on the same procedure in Example 1. Finally, $y_{i}$ is
generated as
$y_{i}=\beta\phi(x_{i})+\bm{\alpha}^{\intercal}\bm{z}_{i}+\epsilon_{i}$. Here
$\bm{\alpha}=(1,1,1,1,1,0,\cdots,0)\in\mathbb{R}^{p}$ indicates that the first
five elements are invalid IVs. We examine four cases: (i) $\beta=0$, (ii)
$\beta=.03$, (iii) $\beta=.05$, (iv) $\beta=.10$. Similarly, (i) is for Type I
error analysis, and (ii)-(iv) are for power analysis, and we construct the
confidence interval for (iii) $\beta=.05$. All empirical results are
summarized in Figure 3 (testing) and Table 2 (CI) based on $n=10000,p=50$, and
$\nu=0.0,0.5$. Moreover, we further consider invalid IVs with correlated
pleiotropy, where $\bm{\theta}=\bm{\theta}_{0}+\bm{\mu}$ and
$\bm{\alpha}=\bm{\alpha}_{0}+\bm{\mu}$ where $\bm{\theta}_{0}$ and
$\bm{\alpha}_{0}$ are simulated with the same procedure in Example 1, and
$\bm{\mu}=(\mu_{1},\cdots,\mu_{5},0,\cdots)^{\intercal}$ with $\mu_{j}\sim
N(0,1)$ for $j=1,\cdots,5$. All empirical results are summarized in Figure A.2
(testing) based on $n=10000,p=50$, and $\nu=0.5$.

Example 3 (Categorical IVs). Note that the proposed method requires that the
instrument variables follow an elliptically symmetric distribution, which is
usually invalid for categorical data. Yet, in practice, a categorical
instrumental variable is often involved in causal inference, say SNP data in
genetic data analysis. In this example, we examine if the proposed method can
be applied to categorical instrumental variables. Specifically, the
instrumental variables $(\bm{z}_{i})_{i=1,\cdots,n}$ are generated as
$\bm{z}_{i}=\bm{\tau}_{i}+\bm{\tau}_{i}^{\prime}$ to mimic the SNP data, where
$\bm{\tau}_{i}$ and $\bm{\tau}_{i}^{\prime}$ are independent Bernoulli trials,
each with a probability of success $0.3$. Moreover, we set $\bm{\theta}\sim
N(\bm{0},\bm{I}_{p})$ and normalize it by its norm, then $x_{i}$ and $y_{i}$
are generated following the same procedure in Example 1. All empirical results
are summarized in Table A.1 (testing), Table A.2 (CI), and Figure A.2
(boxplot) of Appendix A.4.

Example 4 (Weak IVs). In this example, we examine the performance and
stability of the proposed method with weak IVs. Specifically, we set
$\bm{\theta}\sim N(\bm{0},\bm{I}_{p})$, $\theta_{j}=0;j=1,\cdots,\lfloor\pi
p\rfloor$, and normalize it by its norm, then $x_{i}$ and $y_{i}$ are
generated following the same procedure in Example 1 based on $(n=5000,p=50)$,
and $\pi=0.0,0.1,0.3$. All empirical results are summarized in Figure A.3
(testing), Table A.3 (CI) of Appendix A.5.

Example 5 (Non-additive and epistatic effects). In this example, we examine
the performance and stability of the proposed method under non-additive and
epistatic genetic effects. First, $(\bm{z}_{i})_{i=1,\cdots,n}$ are generated
based on the same setting in Example 3. To incorporate the non-additive and
epistatic effects,
$x_{i}=\phi^{-1}(\bm{\theta}^{\intercal}_{a}I(\bm{z}_{i}=1)+\bm{\theta}^{\intercal}_{d}I(\bm{z}_{i}=2)+\sum_{(j,j^{\prime})\in\mathcal{J}}\delta_{j,j^{\prime}}z_{ij}z_{ij^{\prime}}+w_{i})$.
Here, we set $\bm{\theta}_{a}\sim N(\bm{0},\bm{I}_{p})$, and
$\bm{\theta}_{d}=\lambda\bm{\theta}_{a}$ presents non-additive effects when
$\lambda\neq 2$. Besides, $\bm{\delta}\sim
N(\bm{0},0.1\bm{I}_{|\mathcal{J}|})$ presents epistatic (i.e. interaction)
effects, and $\mathcal{J}$ is a set of randomly selected pairs, where each
pair is uniformly sampled. Finally, $y_{i}$ is generated following the same
procedure in Example 3. In this example, we set $n=5000,p=50$,
$\lambda=0.3,0.5$ and $|\mathcal{J}|=\lfloor 0.1p\rfloor,\lfloor 0.3p\rfloor$.
All empirical results are summarized in Figure A.4 (testing), Table A.4 (CI)
of Appendix A.5.

Example 6 (Misspecified models). We examine the performance and stability of
the proposed method for misspecified models. Specifically,
$(\bm{z}_{i},x_{i})_{i=1,\cdots,n}$ are generated with the same procedure in
Example 1. In Stage 2, we consider misspecified models:
$y_{i}=\beta\psi(x_{i})+\epsilon_{i}$ with $\psi(x)=x$, $\psi(x)=\exp(x)$,
$\psi(x)=|x|$, $\psi(x)=1/x$, and $\psi(x)=\log(|x|)$. According to the
simulation results in Example 1, we mainly consider $\phi(x)=x^{2}$ and
$\phi(x)=1/x$ to highlight the differences between the proposed methods and
other competitors. All empirical results are summarized in Figure A.5
(testing) of Appendix A.5.

Overall, the simulation study indicates that the proposed 2SIR method performs
well under all nonlinear transformations to control Type I error, yielding
high power, and constructing valid CIs based on continuous/categorical IVs
with or without invalid IVs. The major conclusions are listed as follows.

  * •

For testing, as suggested in Figures 2 and 3, and A.1, the proposed 2SIR and
its Cauchy-combined test yield competitive performance for “linear”, “cube-
root” and “PL” cases compared with 2SLS and PT-2SLS; and superior performance
for “log”, “inverse”, and “quad” cases. The largest percentage of improvement
in terms of power are 1433% and 820% over 2SLS and PT-2SLS, respectively.

  * •

For CI, as indicated in Tables 1, 2, and A.2, 2SLS and PT-2SLS fail to provide
valid CIs when “inverse” and “quad” transformations are used. For other cases,
the proposed 2SIR yields competitive performance. In general, 2SIR is the only
one that can provide a valid CI under an unknown nonlinear transformation.

  * •

As suggested in Figures 3 and A.3, and Tables 2, A.1 and A.2, the proposed
2SIR continues to perform well with invalid, weak or categorical IVs. As
indicated in Figure A.4 and Table A.4, 2SIR is the most robust method against
dominance and epistatic effects, while PT-2SLS is the most sensitive.

  * •

As indicated in Figure A.5, the proposed methods can control Type I errors,
and are more powerful than the competitors when the transformation is
misspecified.

Figure 2: Empirical Type I error ($\beta_{0}=0$) and power
($\beta_{0}=0.05,0.10,0.15$) of the proposed nonlinear causal test for the
simulated example (marginal effect inference) in Example 1 of Section 3.1.
$(n,p)=(2000,10),(2000,50),(5000,50),(5000,100)$ from up to bottom. Here 2SLS,
PT-2SLS, 2SIR, and Comb-2SIR denote two-stage least square, Yeo-Johnson power
transformed two-stage least square, the proposed method, and the Cauchy
combined proposed method, respectively.

|  | 2SLS |  |  | PT-2SLS |  |  |  2SIR (proposed)
---|---|---|---|---|---|---|---|---
$(n,p)$ |  | coverage | length |  |  | coverage | length |  |  | coverage | length
$(2000,10)$ | linear | 0.944 | 0.132 |  |  | 0.943 | 0.132 |  |  | 0.967 | 0.138
| log | 0.946 | 156.422 |  |  | 0.946 | 0.133 |  |  | 0.925 | 0.136
| cube-root | 1.000 | 0.390 |  |  | 1.000 | 0.436 |  |  | 0.975 | 0.138
| inverse | 0.964 | 0.522 |  |  | 0.930 | 0.134 |  |  | 0.979 | 0.138
| PL | 0.950 | 0.134 |  |  | 0.949 | 0.134 |  |  | 0.971 | 0.138
| quad | 0.831 | 0.093 |  |  | 0.823 | 0.092 |  |  | 0.951 | 0.139
$(2000,50)$ | linear | 0.941 | 0.128 |  |  | 0.943 | 0.129 |  |  | 0.974 | 0.136
| log | 1.000 | 176.916 |  |  | 0.913 | 0.123 |  |  | 0.935 | 0.136
| cube-root | 1.000 | 0.328 |  |  | 0.940 | 0.132 |  |  | 0.976 | 0.136
| inverse | 0.990 | 0.149 |  |  | 0.882 | 0.096 |  |  | 0.979 | 0.131
| PL | 0.943 | 0.126 |  |  | 0.944 | 0.127 |  |  | 0.982 | 0.134
| quad | 0.743 | 0.084 |  |  | 0.743 | 0.083 |  |  | 0.976 | 0.134
$(5000,50)$ | linear | 0.950 | 0.094 |  |  | 0.952 | 0.095 |  |  | 0.978 | 0.095
| log | 1.000 | 95.559 |  |  | 1.000 | 0.090 |  |  | 0.972 | 0.097
| cube-root | 1.000 | 0.215 |  |  | 0.999 | 0.095 |  |  | 0.982 | 0.097
| inverse | 0.801 | 0.209 |  |  | 0.640 | 0.060 |  |  | 0.972 | 0.096
| PL | 0.951 | 0.096 |  |  | 0.960 | 0.096 |  |  | 0.977 | 0.096
| quad | 0.522 | 0.052 |  |  | 0.523 | 0.051 |  |  | 0.976 | 0.095
$(10000,50)$ | linear | 0.952 | 0.074 |  |  | 0.955 | 0.074 |  |  | 0.958 | 0.074
| log | 1.000 | 76.293 |  |  | 0.925 | 0.072 |  |  | 0.955 | 0.074
| cube-root | 0.949 | 0.174 |  |  | 0.960 | 0.074 |  |  | 0.944 | 0.074
| inverse | 0.532 | 0.126 |  |  | 0.380 | 0.042 |  |  | 0.969 | 0.075
| PL | 0.954 | 0.075 |  |  | 0.956 | 0.075 |  |  | 0.959 | 0.075
| quad | 0.287 | 0.040 |  |  | 0.291 | 0.043 |  |  | 0.931 | 0.074

Table 1: Empirical coverage and length of the confidence interval for the
simulated example (marginal effect inference) in Example 1 of Section 3.1.

Figure 3: Empirical Type I error ($\beta_{0}=0$) and power
($\beta=0.05,0.10,0.15$) for the simulated example (invalid IVs) in Example 2
of Section 3.1. $\nu=0.0,0.5$ from top to bottom.

|  | 2SLS |  |  | PT-2SLS |  |  |  2SIR (proposed) |
---|---|---|---|---|---|---|---|---|---
$\nu$ |  | coverage | length |  |  | coverage | length |  |  | coverage | length |
$0.0$ | linear | 0.945 | 0.078 |  |  | 0.945 | 0.078 |  |  | 0.948 | 0.078 |
| log | 0.999 | 79.988 |  |  | 0.928 | 0.078 |  |  | 0.952 | 0.078 |
| cube-root | 0.965 | 0.190 |  |  | 0.972 | 0.079 |  |  | 0.949 | 0.079 |
| inverse | 0.598 | 0.159 |  |  | 0.510 | 0.050 |  |  | 0.954 | 0.078 |
| PL | 0.951 | 0.079 |  |  | 0.950 | 0.079 |  |  | 0.951 | 0.079 |
| quad | 0.443 | 0.043 |  |  | 0.456 | 0.043 |  |  | 0.964 | 0.079 |
$0.5$ | linear | 0.951 | 0.050 |  |  | 0.951 | 0.050 |  |  | 0.945 | 0.050 |
| log | 1.000 | 213.678 |  |  | 0.948 | 0.056 |  |  | 0.946 | 0.050 |
| cube-root | 1.000 | 0.216 |  |  | 0.945 | 0.051 |  |  | 0.943 | 0.049 |
| inverse | 0.827 | 0.210 |  |  | 0.645 | 0.062 |  |  | 0.940 | 0.050 |
| PL | 0.942 | 0.050 |  |  | 0.912 | 0.049 |  |  | 0.936 | 0.050 |
| quad | 0.541 | 0.055 |  |  | 0.514 | 0.055 |  |  | 0.946 | 0.049 |

Table 2: Empirical coverage and length of the confidence interval for the
simulated example (invalid IVs) in Example 2 of Section 3.1.

###  3.2 Simulation for transformation estimation

This subsection examines the proposed adjusted inverse regression (2SIR+AIR)
in (15) under various nonlinear transformations, and the estimation accuracy
is measured by mean square error (MSE) and uniform error (UE):

$\text{MSE}(\widehat{\phi},\phi_{0})=\operatorname{E}\Big{(}\big{(}\widehat{\phi}(x)-\phi_{0}(x)\big{)}^{2}\Big{)},\qquad\text{UE}(\widehat{\phi},\phi_{0})=\operatorname{E}\sup_{x\in\mathcal{X}}\big{|}\widehat{\phi}(x)-\phi_{0}(x)\big{|}$
(15)

where $\mathcal{X}$ is a region of causal interest, which is replaced as 100
grid points of [5%-quantile, 95%-quantile] of $x$ for evaluation. We also
compare the results with a conditional mean function to highlight the role of
the ratio-correction in (6).

Specifically, we set $\bm{\theta}=(p^{-1/2},\cdots,p^{-1/2})^{\intercal}$ and
$\beta=1$ in (1). Note that $\mathcal{D}_{1}$ and $\mathcal{D}_{2}$ are
generated with the same setting in Example 1 in Section 3.1 with
$w_{i}=u_{i}+\gamma_{i}$, $u_{i}$ and $\gamma_{i}$ are independently generated
from $N(0,1)$. Five nonlinear transformations are considered: (1) linear:
$\phi(x)=x$; (2) logarithm: $\phi(x)=\log(x)$; (3) cube root:
$\phi(x)=x^{1/3}$, (4) piecewise linear (PL): $\phi(x)=xI(x\leq
0)+0.5xI(x>0)$, (5) quadratic (quad): $\phi(x)=x^{2}$. Note that the
conditional mean regression (5) is conducted based on a KNN model with the
number of neighbors as 100. The simulation is replicated 100 times with
$n=2000,p=10,50,100$, the resulting MSEs and UEs for all methods are
summarized in Table 3, and the fitted transformations for $p=10$ is
illustrated in Figure 4.

It is evident that the proposed 2SIR+AIR method substantially outperforms 2SLS
and PT-2SLS in most cases, except that 2SLS yields better performance in the
“linear” case where the proposed model in (1) becomes a linear structural
equation model. For other cases, the amount of improvement is significant,
with the largest improvement of (MSE: 99.9%, UE: 96.6%) and (MSE: 92.2%, UE:
64.4%) over 2SLS and PT-2SLS, respectively.

|  | 2SLS |  |  | PT-2SLS |  |  | Cond-mean(KNN) |  |  |  2SIR+AIR (proposed)
---|---|---|---|---|---|---|---|---|---|---|---
$p$ |  | MSE | UE |  |  | MSE | UE |  |  | MSE | UE |  |  | MSE | UE
$10$ | linear | 0.000(.000) | 0.000(.000) |  |  | 0.525(.005) | 1.216(.005) |  |  | 3.530(.179) | 3.076(0.090) |  |  | 0.117(.003) | 0.615(.012)
| log | 363.405(48.756) | 9.892(0.045) |  |  | 0.619(.004) | 1.362(.004) |  |  | 3.471(.205) | 2.945(.094) |  |  | 0.118(.002) | 0.589(.016)
| cube-root | 346.575(6.023) | 21.777(0.042) |  |  | 1.293(.009) | 1.737(.009) |  |  | 3.336(.205) | 2.766(.099) |  |  | 0.113(.002) | 0.584(.016)
| PL | 1.026(.005) | 2.130(.002) |  |  | 0.540(.004) | 1.284(.005) |  |  | 2.853(.207) | 2.614(.104) |  |  | 0.123(.003) | 0.645(.016)
| quad | 2.461(.009) | 3.073(.004) |  |  | 2.083(.009) | 2.824(.004) |  |  | 1.323(.060) | 1.568(.042) |  |  | 0.123(.004) | 0.638(.013)
$50$ | linear | 0.000(.000) | 0.000(.000) |  |  | 0.535(.005) | 1.171(.004) |  |  | 3.305(.214) | 3.022(.096) |  |  | 0.125(.003) | 0.598(.015)
| log | 223.565(22.881) | 12.106(.028) |  |  | 0.616(.004) | 1.342(.003) |  |  | 3.273(.214) | 2.829(.104) |  |  | 0.124(.002) | 0.534(.016)
| cube-root | 355.761(5.317) | 19.961(.038) |  |  | 1.302(.010) | 1.738(.008) |  |  | 3.408(.216) | 2.922(.100) |  |  | 0.121(.003) | 0.561(.013)
| PL | 1.022(.004) | 2.134(.002) |  |  | 0.546(.005) | 1.256(.005) |  |  | 3.113(.214) | 2.965(.100) |  |  | 0.119(.003) | 0.583(.016)
| quad | 2.474(.009) | 3.287(.003) |  |  | 2.095(.008) | 3.033(.004) |  |  | 1.162(.069) | 1.581(.055) |  |  | 0.163(.006) | 0.837(.020)
$100$ | linear | 0.000(.000) | 0.000(.000) |  |  | 0.526(.004) | 1.204(.004) |  |  | 3.203(.217) | 3,019(.095) |  |  | 0.142(.003) | 0.570(.010)
| log | 615.467(32.895) | 7.429(.044) |  |  | 0.623(.005) | 1.580(.004) |  |  | 3.591(.220) | 2.741(.111) |  |  | 0.148(.003) | 0.539(.011)
| cube-root | 354.663(5.198) | 20.571(.023) |  |  | 1.300(.009) | 1.740(.010) |  |  | 3.818(.217) | 3.157(.104) |  |  | 0.140(.003) | 0.565(.012)
| PL | 1.018(.005) | 2.103(.002) |  |  | 0.541(.004) | 1.176(.004) |  |  | 3.638(.219) | 3.057(.107) |  |  | 0.142(.003) | 0.572(.015)
| quad | 2.468(.009) | 3.097(.004) |  |  | 2.092(.008) | 2.851(.004) |  |  | 1.201(.076) | 1.492(.057) |  |  | 0.232(.009) | 1.015(.023)

Table 3: Mean square error (MSE) and uniform error (UE) (standard errors in
parentheses) for the simulated example in Section 3.2. Here cond-mean(KNN),
and 2SIR+AIR denote nonparametric regression in (5), and the proposed method
in (6), respectively.

Figure 4: Fitted transformations of the simulated example in Section 3.2,
where the true transformations are: (1,1) linear; (1,2) logarithm; (2,1) cubic
root; (2,2) quadratic.

##  4 Real data analysis

In this section, we implement the proposed method for an analysis of the
Alzheimer’s Disease Neuroimaging Initiative (ADNI) dataset and the
International Genomics of Alzheimer’s Project (IGAP; [20]) GWAS summary
dataset to identify putative causal genes for Alzheimer’s Disease (AD).
Specifically, the ADNI dataset consists of 819 individual-level subjects (229
cognitively normal, 398 with MCI, and 192 with AD), 17,201 candidate genes,
and 620,901 SNPs, with standard cognitive and functional measures typical of
clinical trials over 12 months. The IGAP GWAS dataset consists of summary
statistics of about 7 million SNPs to AD based on 54,162 samples.

###  4.1 Data preprocessing

To facilitate the analysis, we pre-process the dataset and refine the
candidate SNPs as follows. For the ADNI dataset, we first exclude SNPs with
MAF $\leq$ 0.05, with missing values, or failing the Hardy-Weinberg
equilibrium test at the significant level of 0.001. Next, we further prune the
SNPs to ensure that any of their pairwise Pearson correlations in absolute
values was no more than 0.6. For the IGAP GWAS dataset, we conduct imputation
for summary statistics of missing SNPs by using the software ImpG [32], based
on 489 unrelated individuals with European ancestry from the 1000 Genomes
Project [37], yet remove the imputed SNPs with imputation accuracy smaller
than 0.3. Finally, we define the cis-region of the gene by expanding 100kb
upstream and downstream of its coding region, and take the top 50 intersecting
SNPs (available both on the ADNI dataset and imputed IGAP dataset), with the
largest absolute correlations with the gene’s expression level. Taken
together, the pre-processed dataset consists of 712 individual-level genotype
and gene-expression data with 50 SNPs, and independent summary statistics for
the associated SNPs based on 54,162 samples.

Figure 5: The bar-plot for negative logarithm of $p$-values of significant
genes for AD by at least one method, where the $x$-axis represents genes, the
$y$-axis represents $-\log_{10}(p)$.

###  4.2 Results

Next, 2SLS, PT-2SLS, and the proposed Comb-2SIR are applied to the pre-
processed data. The Quantile-Quantile plots of the negative logarithm p-values
of 17,201 genes for all methods are illustrated in Figure 6. Overall, it is in
agreement with those of other AD studies in the literature [40, 12]. For power
analysis, as indicated in Figure 5, with the Bonferroni adjusted significance
cutoff $0.05/17201$, 20 genes are identified as significantly related to AD by
at least one method. Specifically, among them 12 were significant by 2SLS
and/or PT-2SLS, 18 are significant by SIR and/or Comb-2SIR. Two genes, APOE
and TOMM40 on chromosome 19, are well-known to be related to AD [2, 26, 23];
the former is identified by all three methods while the latter is only
identified by Comb-2SIR. Besides TOMM40, 7 genes, BCL3, BIN1, CEACAM19,
CHRNA2, CLPTM1, FBXO46, NKPD1, are only identified by Comb-2SIR. We searched
these 7 genes in large-scale GWAS results and found all of them except FBXO46
contained genetic variants that have been reported to be significantly
associated with AD [17, 25, 1]. A further literature search gives more
supporting evidence about their associations with AD. Specifically, BCL3 has
been discovered to be associated with late-onset familial AD [29, 34]; in AD
brains, BIN1 has increased expression levels [6, 3]; CEACAM19 has been
suggested as a candidate gene related to human aging [8]; CHRNA2 has been
implicated in potentially contributing to learning and memory functions [30],
and as a potential target of clinical AD drugs [5].

Figure 6: QQ-plots for negative logarithm of p-values of 17,201 genes for AD.

For illustration, Figure 7 shows the fitted transformations for two genes:
APOC1 and BCL3 (others are included in Supplementary), where the former is
identified by all three methods, while the latter is only identified by the
proposed method. As indicated in Figure 7, for APOC1, which is successfully
detected by 2SLS, the estimated transformation by the proposed method is
roughly in agreement with the linear pattern estimated by 2SLS. For BCL3, in
contrast, the estimated transformation by the proposed method is significantly
different from that of 2SLS and PT-2SLS, indicating that the linear pattern
might be invalid here. This may be a reason for less significance given by
2SLS and PT-2SLS, offering practical and empirical evidence for nonlinear
causal effects in a real dataset.

Figure 7: Fitted transformations of two illustrative genes in the real dataset
in Section 4.1. Left panel is for gene APOC1, which is identified by all three
methods. Right panel is for BCL3, which is only identified by the proposed
method.

###  4.3 Negative control outcomes

In this section, we demonstrate Type I error control of 2SLS, PT-2SLS, and the
proposed 2SIR methods based on the ADNI dataset with negative control
outcomes. Specifically, we implement the methods based on individual-level
SNPs and gene expressions while generating negative control outcomes by
simulating random noises from a standard normal distribution so that no gene
is causal to the outcome. In this case, the p-values in the QQ-plot are
expected to follow a uniform distribution. Figure 8 exhibits the QQ-plots of
$-\log_{10}(p)$ of the methods, suggesting that the p-values provided by 2SLS,
PT-2SLS, and 2SIR are appropriately distributed in this negative control
dataset.

Figure 8: QQ-plots for 2SLS, PT-2SLS, and the proposed 2SIR methods based on
ADNI negative control outcomes.

##  5 Discussion and conclusions

Nonlinear modeling in TWAS has broad applicability and potential significance
in identifying causal exposure-outcome associations, such as gene-trait
associations. However, nonlinear causal inference is plagued by the lack of
individual-level GWAS data (with only summary statistics for the outcome
available). In this article, we have proposed a flexible causal model for GWAS
summary data while allowing an arbitrary nonlinear causal effect of exposure,
substantially relaxing the common assumption of linearity in the current
practice of TWAS. A novel method called 2SIR+AIR is developed to estimate the
marginal causal effect and the nonlinear transformation. In addition, we have
developed inferential tools to assess exposure-outcome associations, including
hypothesis testing and interval estimation.

We have demonstrated the applicability of the proposed model and methods by
studying the ADNI gene expression and the IGAP AD GWAS datasets to identify
putative causal genes for AD. Our results suggest that the proposed method
agrees with two existing methods (2SLS and PT-2SLS) in 10 of 12 putative
causal genes, but it additionally identifies 8 other genes (TOMM40, BCL3,
BIN1, CEACAM19, CHRNA2, CLPTM1, FBXO46, NKPD1). We do find evidence supporting
associations between these additional genes and AD, especially TOMM40 (a well-
known AD-related gene) We also observe higher $R^{2}$’s for the stage one
model of the proposed method than the $R^{2}$’s based on existing linear
models, offering another source of evidence that nonlinear causal effects are
likely to be present in real data. Furthermore, the differences in the
estimated nonlinear transformations among 8 potential genes (with their
figures included in Supplementary) from the estimated linear patterns by 2SLS
are substantial, and not explained by previous works. Our finding reasonably
suggests potential nonlinearity in gene-trait causal associations based on
GWAS data. We believe that the proposed method has great potential and could
further advance research in TWAS, including nonlinear treatment effect
analysis, subgroup analysis, and robustness analysis. Finally, in addition to
TWAS, the proposed method can be equally applied to study other exposure-
outcome causal relationships in a more general context.

## References

  * [1] G. W. Beecham, K. Hamilton, A. C. Naj, E. R. Martin, M. Huentelman, A. J. Myers, J. J. Corneveaux, J. Hardy, J.-P. Vonsattel, S. G. Younkin, et al.  Genome-wide association meta-analysis of neuropathologic features of Alzheimer’s disease and related dementias.  PLoS Genetics, 10(9):e1004606, 2014.
  * [2] G. Bu.  Apolipoprotein E and its receptors in Alzheimer’s disease: pathways, pathogenesis and therapy.  Nature Reviews Neuroscience, 10(5):333–344, 2009.
  * [3] J. Chapuis, F. Hansmannel, M. Gistelinck, A. Mounier, C. Van Cauwenberghe, K. Kolen, F. Geller, Y. Sottejeau, D. Harold, P. Dourlen, et al.  Increased expression of BIN1 mediates Alzheimer genetic risk by modulating tau pathology.  Molecular Psychiatry, 18(11):1225–1234, 2013.
  * [4] R. D. Cook.  Regression Graphics: Ideas for Studying Regressions Through Graphics, volume 482.  John Wiley & Sons, 2009.
  * [5] J. Cummings, G. Lee, A. Ritter, M. Sabbagh, and K. Zhong.  Alzheimer’s disease drug development pipeline: 2019.  Alzheimer’s & Dementia: Translational Research & Clinical Interventions, 5:272–293, 2019.
  * [6] P. L. De Jager, G. Srivastava, K. Lunnon, J. Burgess, L. C. Schalkwyk, L. Yu, M. L. Eaton, B. T. Keenan, J. Ernst, C. McCabe, et al.  Alzheimer’s disease: early alterations in brain DNA methylation at ANK1, BIN1, RHBDF2 and other loci.  Nature Neuroscience, 17(9):1156–1163, 2014.
  * [7] N. Duan and K.-C. Li.  Slicing regression: a link-free regression method.  Annals of Statistics, 19(2):505–530, 1991.
  * [8] D. S. Evans and S. R. Cummings.  Identification of ADAMTS7 and CEACAM19 as candidate healthy aging associated genes.  Innovation in Aging, 3(Supplement_1):S102–S102, 2019.
  * [9] J. Fan and R. Li.  Variable selection via nonconcave penalized likelihood and its oracle properties.  Journal of the American Statistical Association, 96(456):1348–1360, 2001.
  * [10] J. H. Friedman.  Greedy function approximation: a gradient boosting machine.  Annals of Statistics, 29(5):1189–1232, 2001.
  * [11] E. R. Gamazon, H. E. Wheeler, K. P. Shah, S. V. Mozaffari, K. Aquino-Michaels, R. J. Carroll, A. E. Eyler, J. C. Denny, D. L. Nicolae, N. J. Cox, et al.  A gene-based association method for mapping traits using reference transcriptome data.  Nature Genetics, 47(9):1091–1098, 2015.
  * [12] J. Gockley, K. S. Montgomery, W. L. Poehlman, J. C. Wiley, Y. Liu, E. Gerasimov, A. K. Greenwood, S. K. Sieberts, A. P. Wingo, T. S. Wingo, et al.  Multi-tissue neocortical transcriptome-wide association study implicates 8 genes across 6 genomic loci in Alzheimer’s disease.  Genome Medicine, 13(1):1–15, 2021.
  * [13] Z. Guo, H. Kang, T. Tony Cai, and D. S. Small.  Confidence intervals for causal effects with invalid instruments by using two-stage hard thresholding with voting.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 80(4):793–815, 2018.
  * [14] A. Gusev, A. Ko, H. Shi, G. Bhatia, W. Chung, B. W. Penninx, R. Jansen, E. J. De Geus, D. I. Boomsma, F. A. Wright, et al.  Integrative approaches for large-scale transcriptome-wide association studies.  Nature Genetics, 48(3):245–252, 2016.
  * [15] P. Hall, J. Horowitz, et al.  A simple bootstrap method for constructing nonparametric confidence bands for functions.  Annals of Statistics, 41(4):1892–1921, 2013.
  * [16] A. Inoue and G. Solon.  Two-sample instrumental variables estimators.  The Review of Economics and Statistics, 92(3):557–561, 2010.
  * [17] I. E. Jansen, J. E. Savage, K. Watanabe, J. Bryois, D. M. Williams, S. Steinberg, J. Sealock, I. K. Karlsson, S. Hägg, L. Athanasiu, et al.  Genome-wide meta-analysis identifies new loci and functional pathways influencing Alzheimer’s disease risk.  Nature Genetics, 51(3):404–413, 2019.
  * [18] H. Kang, A. Zhang, T. T. Cai, and D. S. Small.  Instrumental variables estimation with some invalid instruments and its application to Mendelian randomization.  Journal of the American statistical Association, 111(513):132–144, 2016.
  * [19] H. Kang, A. Zhang, T. T. Cai, and D. S. Small.  Instrumental variables estimation with some invalid instruments and its application to Mendelian randomization.  Journal of the American statistical Association, 111(513):132–144, 2016.
  * [20] J.-C. Lambert, C. A. Ibrahim-Verbaas, D. Harold, A. C. Naj, R. Sims, C. Bellenguez, G. Jun, A. L. DeStefano, J. C. Bis, G. W. Beecham, et al.  Meta-analysis of 74,046 individuals identifies 11 new susceptibility loci for Alzheimer’s disease.  Nature Genetics, 45(12):1452–1458, 2013.
  * [21] K.-C. Li.  Sliced inverse regression for dimension reduction.  Journal of the American Statistical Association, 86(414):316–327, 1991.
  * [22] Y. Liu and J. Xie.  Cauchy combination test: a powerful test with analytic p-value calculation under arbitrary dependency structures.  Journal of the American Statistical Association, 115(529):393–402, 2020.
  * [23] D. M. Lyall, S. E. Harris, M. E. Bastin, S. M. Maniega, C. Murray, M. W. Lutz, A. M. Saunders, A. D. Roses, M. d. C. V. Hernández, N. A. Royle, et al.  Alzheimer’s disease susceptibility genes APOE and TOMM40, and brain white matter integrity in the Lothian Birth Cohort 1936.  Neurobiology of Aging, 35(6):1513–e25, 2014.
  * [24] T. F. Mackay.  Epistasis and quantitative traits: using model organisms to study gene–gene interactions.  Nature Reviews Genetics, 15(1):22–33, 2014.
  * [25] R. E. Marioni, S. E. Harris, Q. Zhang, A. F. McRae, S. P. Hagenaars, W. D. Hill, G. Davies, C. W. Ritchie, C. R. Gale, J. M. Starr, et al.  GWAS on family history of Alzheimer’s disease.  Translational Psychiatry, 8(1):1–7, 2018.
  * [26] A. Mise, Y. Yoshino, K. Yamazaki, Y. Ozaki, T. Sao, T. Yoshida, T. Mori, Y. Mori, S. Ochi, J.-i. Iga, et al.  TOMM40 and APOE gene expression and cognitive decline in Japanese Alzheimer’s disease subjects.  Journal of Alzheimer’s Disease, 60(3):1107–1117, 2017.
  * [27] J. Morrison, N. Knoblauch, J. H. Marcus, M. Stephens, and X. He.  Mendelian randomization accounting for correlated and uncorrelated pleiotropic effects using genome-wide summary statistics.  Nature Genetics, 52(7):740–747, 2020.
  * [28] E. A. Nadaraya.  On estimating regression.  Theory of Probability & Its Applications, 9(1):141–142, 1964.
  * [29] K. Nho, S. Kim, E. Horgusluoglu, S. L. Risacher, L. Shen, D. Kim, S. Lee, T. Foroud, L. M. Shaw, J. Q. Trojanowski, et al.  Association analysis of rare variants near the APOE region with CSF and neuroimaging biomarkers of Alzheimer’s disease.  BMC Medical Genomics, 10(1):45–52, 2017.
  * [30] H. Nichol.  Optogenetic Investigation of Chrna2 Cells in The Subiculum and Their Role in Modulating Entorhinal Cortex Input.  McGill University (Canada), 2015.
  * [31] P. C. Okoro, R. Schubert, X. Guo, W. C. Johnson, J. I. Rotter, I. Hoeschele, Y. Liu, H. K. Im, A. Luke, L. R. Dugas, et al.  Transcriptome prediction performance across machine learning models and diverse ancestries.  Human Genetics and Genomics Advances, 2(2):100019, 2021.
  * [32] B. Pasaniuc, N. Zaitlen, H. Shi, G. Bhatia, A. Gusev, J. Pickrell, J. Hirschhorn, D. P. Strachan, N. Patterson, and A. L. Price.  Fast and accurate imputation of summary statistics enhances evidence of functional enrichment.  Bioinformatics, 30(20):2906–2914, 2014.
  * [33] J. Pearl.  Causality.  Cambridge University Press, 2009.
  * [34] M. Pericak-Vance, J. Bebout, P. Gaskell, L. Yamaoka, W.-Y. Hung, M. Alberts, A. Walker, R. Bartlett, C. Haynes, K. Welsh, et al.  Linkage studies in familial Alzheimer disease: evidence for chromosome 19 linkage.  American Journal of Human Genetics, 48(6):1034, 1991.
  * [35] X. Shen and H.-C. Huang.  Grouping pursuit through a regularization solution surface.  Journal of the American Statistical Association, 105(490):727–739, 2010.
  * [36] N. Solovieff, C. Cotsapas, P. H. Lee, S. M. Purcell, and J. W. Smoller.  Pleiotropy in complex traits: challenges and strategies.  Nature Reviews Genetics, 14(7):483–495, 2013.
  * [37] The 1000 Genomes Project Consortium.  A global reference for human genetic variation.  Nature, 526(7571):68, 2015.
  * [38] G. Wahba.  Spline Models for Observational Data.  SIAM, 1990.
  * [39] F. Windmeijer, H. Farbmacher, N. Davies, and G. Davey Smith.  On the use of the lasso for instrumental variables estimation with some invalid instruments.  Journal of the American Statistical Association, 114(527):1339–1350, 2019.
  * [40] A. P. Wingo, Y. Liu, E. S. Gerasimov, J. Gockley, B. A. Logsdon, D. M. Duong, E. B. Dammer, C. Robins, T. G. Beach, E. M. Reiman, et al.  Integrating human brain proteomes with genome-wide association data implicates new proteins in Alzheimer’s disease pathogenesis.  Nature Genetics, 53(2):143–146, 2021.
  * [41] H. Xue, X. Shen, and W. Pan.  Constrained maximum likelihood-based mendelian randomization robust to both correlated and uncorrelated pleiotropic effects.  The American Journal of Human Genetics, 108(7):1251–1269, 2021\.
  * [42] J. Yang, B. Benyamin, B. P. McEvoy, S. Gordon, A. K. Henders, D. R. Nyholt, P. A. Madden, A. C. Heath, N. G. Martin, G. W. Montgomery, et al.  Common SNPs explain a large proportion of the heritability for human height.  Nature Genetics, 42(7):565–569, 2010.
  * [43] I.-K. Yeo and R. A. Johnson.  A new family of power transformations to improve normality or symmetry.  Biometrika, 87(4):954–959, 2000.
  * [44] C.-H. Zhang.  Nearly unbiased variable selection under minimax concave penalty.  Annals of Statistics, 38(2):894–942, 2010.
  * [45] W. Zhang and D. Ghosh.  On the use of kernel machines for Mendelian randomization.  Quantitative Biology, 5(4):368–379, 2017.
  * [46] L.-X. Zhu and K. W. Ng.  Asymptotics of sliced inverse regression.  Statistica Sinica, 5(2):727–736, 1995.
  * [47] Z. Zhu, F. Zhang, H. Hu, A. Bakshi, M. R. Robinson, J. E. Powell, G. W. Montgomery, M. E. Goddard, N. R. Wray, P. M. Visscher, et al.  Integration of summary data from GWAS and eQTL studies predicts complex trait gene targets.  Nature Genetics, 48(5):481–487, 2016.


# Real and virtual direct photon measurements with ALICE

H. Sebastian Scheid
for the ALICE Collaboration
<EMAIL_ADDRESS>

###### Abstract

In this contribution the latest measurements of real and virtual photons in
Pb–Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV from the ALICE Collaboration
are presented. The extracted photon spectra are compared to predictions from
state-of-the-art models that explicitly include the whole evolution of the
heavy-ion collision as well as measurements at lower energies. Experimental
developments that in the future will help to extend the range and procession
of the measurements are introduced before concluding with an outlook to
possible future measurements at the LHC.

##  1 Why Photons?

If we think about the creation of the little big bang in the laboratory by
colliding heavy nuclei at ultra-relativistic velocities we can identify
different stages of these collisions in which different physics is at work.
While most partons in the nuclei will just pass through each other and deposit
energy in the collision zone, some will interact in most violent initial
collisions. These initial scatterings can be described in perturbative
quantum-chromo dynamics (pQCD). After this the system will start to build up a
medium which will approach equilibrium and form a plasma of quarks and gluons
at a thermal scale of some hundred MeV. This can be described by modelling a
pre-equilibrium phase followed by hydrodynamic evolution. With the passing of
time, this system will expand and cool down. Confinement will force the quarks
and gluons into hadrons. Shortly after, the inelastic interactions will cease,
followed by elastic scatterings. Finally the debris of the collision will be
measured in the detector. While hadronic observables will always be obscured
by interactions in the later stages of the collision, electromagnetic probes
carry their information about the early stages of the heavy-ion collision to
the detector unscathed. In addition, a time ordering of the transverse
momentum ($p_{\rm T}$), and in the dielectron case also invariant mass
($m_{\rm ee}$), is expected. Direct photons, that are produced earlier in the
collision, tend to have higher $p_{\rm T}$ or $m_{\rm ee}$. In the picture of
an equilibrated thermal source producing the photons, the spectrum should
depend on the size and the temperature of the source. While the absolute yield
depends on the space time volume of the source, the inverse slope of the
spectrum is only sensitive to the temperature $T$. In the case of real photons
the energy of the photon measured in the detector is blue shifted due to the
expanding source. For the virtual direct photons such a blue shift is not
present since the Lorentz invariant mass of the virtual photons is the
observable. In the terminology this is reflected by measuring an effective
temperature in the real direct photon measurements.

While the above arguments are very compelling, the measurement of real and
virtual direct photons are notoriously challenging. In both measurements the
main sources of real photons and dielectrons are decays of scalar and pseudo-
scalar mesons. In the case of virtual photons an additional physical
background stems from the semi-leptonic decays of correlated open heavy-
flavour (HF) hadrons. Since not all decay products but only the leptons are
reconstructed, the mass distribution of the HF contributions is wide and with
little structure. In addition it is not a priori clear which of the electrons
originate from the same hadronic decay, which induces a combinatorial
background of fake pairs in the dielectron analysis. This combinatorial
background of fake pairs is estimated by pairing $\rm{e^{+}e^{+}}$ and
$\rm{e^{-}e^{-}}$ which are then subtracted from all $\rm{e^{+}e^{-}}$ pairs.

##  2 Results

Discussed here are measurements of photons in the ALICE detector using the
internal and external conversion, by identifying the electron candidates. In
addition ALICE has the possibility to directly measure photons in the
electromagnetic calorimeters, e.g. the photon spectrometer (PHOS). The
presented data from Pb–Pb collisions at a center-of-mass energy per nucleon
pair $\sqrt{s}_{\rm{NN}}=5.02$ TeV was recorded with a minimum bias trigger
(2015), and an online enhancement of events in the 0–10% and 30–50% centrality
classes (2018).

To extract the direct-photon signal in the dielectron measurement a
parameterisation of the invariant mass $m_{\rm ee}$ spectrum with a three
component function is performed: $f_{\rm fit}=r\times f_{\rm dir}+(1-r)\times
f_{\rm LF}+f_{\rm HF}$. Here $f$ denotes the contributions from direct photons
(dir), as well as from the light-flavour (LF) and the heavy-flavour (HF)
hadron decays. The free parameter $r$ reflects, the fraction of virtual direct
to inclusive virtual photons which is equivalent to the same fraction for real
photons for $m_{\rm ee}\rightarrow 0$. The HF contribution is fixed to the
expected contribution. The range of the parametrisation is chosen to be always
above the mass of the $\pi^{0}$-meson. To calculate the yield of real photons
from $r$, the inclusive photon $\gamma^{\rm inc}$ measurement is needed. This
was performed with the photon conversion method (PCM).

Figure 1: Invariant yield (left) of direct photons and $R_{\gamma}$ (right) at
different centralities in Pb–Pb collisions at 5.02 TeV. The data are compared
to pQCD [11] and model calculations that include all stages of the collision
[7].

The measurement of direct photons is shown for different centrality intervals
in fig. 1 (left). The direct photons were measured via reconstructing real
photons in the PCM method and with the virtual direct photon method in the
0–10% centrality class. Upper limits as a 95% confidence interval (CL) are
given in case of data points being consistent with 0 within statistical
uncertainties. The data are compared to model calculations for direct photons
from pQCD processes. In all centrality classes, the data tend to be above the
pQCD prediction at an intermediate $p_{\rm T}$ range of 1-3 GeV/$c$ . A good
agreement of the data in all measured centrality classes is achieved by model
calculations that explicitly include all the different stages of a heavy-ion
collision that emit photons, i.e. a pre-equilibrium as well as a QGP and hot
hadronic phase [7]. The representation of the measurements in terms of the
ratio of inclusive photons to decay photons $R_{\gamma}=\frac{\gamma^{\rm
inc}}{\gamma^{\rm decay}}$ is shown in fig. 1 (right), for the different
centrality classes. The $R_{\gamma}$ shows that at small $p_{\rm T}$ the
signal of direct photons is expected to be in the order of 5% wrt to the decay
contribution for the given system. The measured data are in agreement with
unity, but also the expected photon signal from model calculations [7]. At
larger $p_{\rm T}$ the data differs from unity indicating a significant source
of photons that are not originating from hadronic decays. This signal is most
pronounced in the most central collisions, and fades away when decreasing the
centrality. This trend is well described by the model calculations.

A comparison of the ALICE measurements with previous measurements of direct
photons by the STAR and the PHENIX collaborations is shown in fig. 2. To
compare the data at different collision energies and centralities, the ratios
to the same model [7] are compared. While the data from the virtual and real
photons from ALICE and the real photon measurement by the STAR collaboration
are in agreement with each other and the model calculations, the two
measurements by the PHENIX collaboration are showing a slight tension wrt the
model calculations. However, taking the systematic uncertainty that dominate
both measurements serious, none of the measurements shows a tension with this
state-of-the-art model.

Figure 2: Measured direct photon yields from the STAR, the PHENIX, and ALICE
collaborations are compared in terms of their ratio to a state-of-the-art
model at their respective energies (left). Correlation function $C$ of two
photons as a function of $\rm Q_{inv}$ in the range of $0.25<k_{\rm T}<0.35$
GeV/$c$ in data and a HIJING simulation in 0–10% Pb–Pb collisions (right).

An alternative approach to extract the direct-photon signal is through
Hanbury-Brown-Twiss (HBT) correlations. Following the WA98 collaboration [4]
the correlation function $C$ of two photons from the same source as a function
of the photon pair relative momentum $Q_{inv}$ is calculated as: ${\rm
C(Q_{inv})=Q_{inv,SE}/Q_{inv,ME}}$. Where SE and ME denote pairs from the same
and different (mixed) events, respectively. Figure 2 (right) shows the ALICE
result in the two photon average momentum interval $0.25<k_{\rm T}<0.35$
GeV/$c$ for data and HIJING simulation in 0–10% most central Pb–Pb collisions.
In this analysis one photon was reconstructed with the PCM while the second
was measured with the PHOS. The sharp peak in simulation and data is expected
from the $\pi^{0}\rightarrow\gamma\gamma$ decay. In the case of photons, ie
vanishing mass, $\rm Q_{inv}$ reduces to the invariant mass of the photon pair
$m_{\gamma\gamma}$. A residual background in the data from particle showers is
parametrised with a second order polynomial. At small $\rm Q_{inv}$ a hint for
an additional signal is quantified by fitting the data with eq. 1, including
the correlation strength $\rm\lambda_{inv}$ and the source size $\rm R_{inv}$.

$\rm C(Q_{inv})=1+\lambda_{inv}\exp(-R^{2}_{inv}Q_{inv}^{2})$ (1)

The large statistical and systematic uncertainties unfortunately do not allow
for a conclusion, but the extraction of the correlation function could in the
future be a good tool to distinguish early and late emission of photons [8].

One of the motivations given in the beginning of this article is the
measurement of the average temperature in the early phase of the collision by
extracting the inverse slope of the virtual direct-photon invariant mass
spectrum at masses above the $\phi$ meson. Figure 3 shows the measured
dielectron yield in the ALICE acceptance as a function of the invariant mass.
The top panel shows the measurement compared with the expectation from
hadronic decays, the so-called hadronic cocktail. Two different approaches for
the contribution from heavy-flavour hadron decays are presented. In the first
case the measured cross section of dielectrons from heavy-flavour hadron
decays measured in pp collisions at $\sqrt{s}=5.02$ TeV [2] is scaled with the
number of binary nucleon-nucleon collisions ($N_{\rm coll}$) (solid lines).
The functional shape as a function of mass is based on the event generator
POWHEG. Measurements of the nuclear modification factor $R_{\rm AA}$ of
electrons from charm and beauty quarks however show a significant
modification, of the $p_{\rm T}$ spectra [3]. To successfully describe the
observed modification models need to include initial cold nuclear matter (CNM)
effects as well as later stage energy loss.

Figure 3: Left: Dielectron yield as a function of $m_{\rm ee}$ measured in
0–10% Pb–Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV (top panel). The data
is compared to hadronic cocktail simulations. In one case (solid lines) the
heavy-flavour contribution is based on the measured cross section in pp
collisions [1] and scaled with $N_{\rm coll}$. In the second case (dashed
lines) the measured $R_{\rm AA}$ of electrons from heavy-flavour hadron decays
is included as described in the text. The respective ratios of the data to the
cocktails are shown in the middle and bottom panel. In both cases also
predictions for thermal dielectrons from a hadronic and partonic phase are
indicated [10, 9]. Right: Dielectron yield as a function of $\rm DCA_{ee}$ in
the intermediate mass region. If the data are consistent with zero upper
limits are given. The spectrum was parametrised with templates for charm,
beauty and a prompt contribution.

The problem in propagating this measured modification to the dielectron
spectra is that the modifications from CNM and energy loss effects affect the
dielectron spectra in different ways. While the energy loss affects each of
the quarks, and thus the electrons, independently, the CNM effects affect the
pair. To accommodate for this an approach was chosen that is able to include
both effects with the expected impact on the pair. As a starting point the
measured $R_{\rm AA}$ of $\rm c,b\rightarrow e$ was parameterised. To
disentangle CNM and energy loss (EL) effects a calculation based on the EPS09
nPDF [6] set was used. By dividing the data parameterisation by the nPDF
expectation, one can construct the expectation for the $R_{\rm AA}$ that
originates only from the energy loss of the heavy quark. In a next step two
different $R_{\rm AA}$ for the dielectrons, one for energy loss and one for
CNM effects, can be calculated in a Monte-Carlo approach. Here the single
electron $R_{\rm AA}$ are used as weights. In the energy loss case
factorisation is assumed and the weight for the pair is calculated as the
product of the electron weights $R^{\rm HFee}_{\rm AA,EL}=w_{1}\times w_{2}$.
In the case of CNM effects no factorisation is expected and the weight for the
pair is the average of the electron weights $R^{\rm HFee}_{\rm
AA,CNM}=(w_{1}+w_{2})/2$. The heavy-flavour dielectron nuclear modification
factor can then be calculated as $R^{\rm HFee}_{\rm AA}=R^{\rm HFee}_{\rm
AA,CNM}\times R^{\rm HFee}_{\rm AA,EL}$. This factor is then used to modify
the cocktail based on $N_{\rm coll}$-scaling. The effect can be seen when
comparing the dashed lines in the upper left panel of fig. 3 with the solid
lines. The middle and bottom panels show the ratio of the measurement to the
cocktail expectation for the case of $N_{\rm coll}$ scaling and including the
$R^{\rm HFee}_{\rm AA}$, respectively. In the mass region between the $\phi$
and the $J/\psi$ we can see that the data is below the expectation from
$N_{\rm coll}$-scaling (middle panel) which is expected, since there is an
overall suppression seen in the $R^{\rm c,b\rightarrow e}_{\rm AA}$. Including
the $R^{\rm HFee}_{\rm AA}$ estimate slightly brings down the cocktail which
leads to the cocktail being closer to the measurement (bottom panel). However,
as discussed before, the uncertainties are large, and no conclusion is
possible. Calculations for a dielectron contribution from a hadronic and a
partonic phase [9, 10] are shown in both ratios to indicate the size of the
expected signal.

A different and very promising approach is to separate the dielectrons from
heavy-flavour hadron decays from the prompt thermal contribution based on
their decay topology. The basic idea is simple: the finite decay length of
charm and beauty hadrons can be exploited to separate the contributions of
heavy-flavour and prompt dielectrons. For this the pair distance-of-closest-
approach ($\rm DCA_{ee}$) is constructed as:

$\rm DCA_{ee}=\sqrt{\frac{DCA_{1}^{2}+DCA_{2}^{2}}{2}}.$ (2)

$\rm DCA_{1}$ and $\rm DCA_{2}$ denote the DCA of the tracks to the
reconstructed primary vertex, normalised to the respective resolution. By
construction we expect a ordering of the distributions depending on the decay
length ($c\tau$) of the decaying particle: $\rm
DCA_{ee}(prompt)<DCA_{ee}(charm)<DCA_{ee}(beauty)$ [1]. A big advantage of
this approach is its model agnostics. Figure 3 (right) shows the data in the
mass interval of 1.2 - 2.6 GeV/$c^{2}$. The data was parametrised with
templates for charm, beauty, $J/\psi$, and a prompt contribution that can be
interpreted as a thermal source. While the uncertainties are large, the data
are well described when including the thermal contribution dominant at low
$\rm DCA_{ee}$, and suppression of the heavy-flavour contributions.

##  3 Conclusion and outlook

The measurement of the direct-photon yield in Pb–Pb collisions at
$\sqrt{s_{\rm NN}}=5.02$ TeV in several centrality classes by the ALICE
Collaboration are well described by model calculations when including all
stages of the heavy-ion collision. The comparison of measurements to those
from the STAR and PHENIX collaborations with respect to the model predictions
at the respective energies shows a good agreement of the ALICE and STAR
measurements with the model. At low $p_{\rm T}$ the central values of the
PHENIX results tend to be above the model expectation, however with large
uncertainties. An alternative approach to access the direct-photon
contribution is the measurement of the correlation function of the direct
photons. A proof of concept was presented, but better statistical precision
will be needed to conclude.

To access the average temperature in the earlier stages of the heavy-ion
collision, i.e. the QGP, it is necessary to understand the dielectron
contribution from charmed hadron decays including its modification. An
approach to include these was discussed, however suffers from large
uncertainties due to its model dependences. A second, model agnostic, approach
uses the decay typologies of the respective contributions. Based on the decay
length of the mother particle the charm and beauty contributions can be
singled out and a possible thermal source can be isolated. Both approaches
have their advantages and disadvantages, while the modeling in the cocktail
introduces model dependencies and large uncertainties, the topological
separation however requires a large amount of statistics.

The recent upgrades of the ALICE TPC and ITS grant a better pointing
resolution (factor 3-6) and increased readout rate (factor 100). Both greatly
benefit the measurements of real and virtual photons. For the 2030s the ALICE
collaboration foresees the construction of a next-generation heavy-ion
experiment [5]. The ALICE 3 detector concept is based around a high-acceptance
all-silicone tracker with unprecedented precision, complemented by particle
identification. Its possibility to efficiently reject HF contributions to the
dielectron spectrum together with the boost in effective statistics will
enable multi-differential analyses. These types of measurements are expected
to be sensitive to the time evolution of the QGP, and its early phases.

## References

  * [1] Shreyasi Acharya et al.  Dielectron production in proton-proton collisions at $\sqrt{s}=7$ TeV.  JHEP, 09:064, 2018.
  * [2] Shreyasi Acharya et al.  Dielectron production in proton-proton and proton-lead collisions at $\sqrt{s_{NN}}=$ 5.02 TeV.  Phys. Rev. C, 102(5):055204, 2020.
  * [3] Shreyasi Acharya et al.  Measurement of electrons from semileptonic heavy-flavour hadron decays at midrapidity in pp and Pb-Pb collisions at $\sqrt{s_{\rm{NN}}}$ = 5.02 TeV.  Phys. Lett. B, 804:135377, 2020.
  * [4] M. M. Aggarwal et al.  Interferometry of direct photons in central Pb-208+Pb-208 collisions at 158-A-GeV.  Phys. Rev. Lett., 93:022301, 2004.
  * [5] ALICE.  Letter of intent for ALICE 3: A next generation heavy-ion experiment at the LHC.  Technical report, CERN, Geneva, 2022.
  * [6] K. J. Eskola, H. Paukkunen, and C. A. Salgado.  EPS09: A New Generation of NLO and LO Nuclear Parton Distribution Functions.  JHEP, 04:065, 2009.
  * [7] Charles Gale, Jean-François Paquet, Björn Schenke, and Chun Shen.  Multimessenger heavy-ion collision physics.  Phys. Rev. C, 105(1):014909, 2022.
  * [8] Oscar Garcia-Montero, Nicole Löher, Aleksas Mazeliauskas, Jürgen Berges, and Klaus Reygers.  Probing the evolution of heavy-ion collisions using direct photon interferometry.  Phys. Rev. C, 102(2):024915, 2020.
  * [9] Ralf Rapp.  Dilepton Spectroscopy of QCD Matter at Collider Energies.  Adv. High Energy Phys., 2013:148253, 2013.
  * [10] Taesoo Song, Wolfgang Cassing, Pierre Moreau, and Elena Bratkovskaya.  Open charm and dileptons from relativistic heavy-ion collisions.  Phys. Rev. C, 97(6):064907, 2018.
  * [11] W. Vogelsang and M. R. Whalley.  A Compilation of data on single and double prompt photon production in hadron hadron interactions.  J. Phys. G, 23:A1–A69, 1997.