# Bank transactions embeddings help to uncover current macroeconomics

1st Maria Begicheva
CDISE
Skolkovo Institute of Science and Technology
Moscow, Russia
<EMAIL_ADDRESS>2nd Alexey Zaytsev
CDISE
Skolkovo Institute of Science and Technology
Moscow, Russia
<EMAIL_ADDRESS>
###### Abstract

Macroeconomic indexes are of high importance for banks: many risk-control
decisions utilize these indexes. A typical workflow of these indexes
evaluation is costly and protracted, with a lag between the actual date and
available index being a couple of months. Banks predict such indexes now using
autoregressive models to make decisions in a rapidly changing environment.
However, autoregressive models fail in complex scenarios related to
appearances of crises.

We propose to use clients’ financial transactions data from a large Russian
bank to get such indexes. Financial transactions are long, and a number of
clients is huge, so we develop an efficient approach that allows fast and
accurate estimation of macroeconomic indexes based on a stream of transactions
consisting of millions of transactions. The approach uses a neural networks
paradigm and a smart sampling scheme.

The results show that our neural network approach outperforms the baseline
method on hand-crafted features based on transactions. Calculated embeddings
show the correlation between the client’s transaction activity and bank
macroeconomic indexes over time.

###### Index Terms:

Bank transactions, macroeconomics, deep learning

##  I Introduction

The coronavirus pandemic has caused great damage to the health and changed
lives of many millions of people around the world. As the state tried its best
to stop the further spread of the virus or at least slow it down, their
actions caused economic problems. Many industries suffered, affecting the
banking industry. In particular, machine learning models degraded due to data
distribution shift and other reasons [1]. Models for assessing the loan
portfolio also suffered, as they always do during crises [2, 3]. The
shortcomings of these models surfaced and the question arose whether it was
worth trusting them at all, if they provide little help during vital crisis
times for banks.

Existing credit risk models in bank are of two types: one type describes
client behaviour at the application stage and following lifetime of this
particular credit, second type describes lifetime of the whole portfolio.
Models of the first type use the loan application form and the client’s credit
history as data. These models use machine learning techniques such as logistic
regression or decision trees to estimate a customer’s probability of default
(PD), which is the probability of repaying a debt to the bank. Moreover, novel
models use deep learning methods, for example, Embedding-Transactional
Recurrent Neural Network (E.T.-RNN) rely on history of bank customers’ credit
and debit card transactions to compute PDs [4].

Both types of models suffer from the change of macroeconomics. For e.g. some
borrower is identified as a ”good” one by model, but when COVID-19 happened
this borrower defaulted on a debt, because he lost his job and was unable to
pay his payments. Second type of models are linear multivariate autoregressive
models, that forecast several different important indexes for bank. There is
no correction for macroeconomics in these models also.

In the current work we consider an important bank’s macroeconomic index: a
default rate. The default rate is a fraction of bank clients, who went overdue
90+ days for the current date [5]. The main goal of this work is to show the
correlation between people transaction activity and default rate over time.
The scheme of our approach is presented in Figure 1.

Our main claims are the following:

  * •

We proposed an approach that can work with the data of complex structure: we
have a time series, and each object in this time series is a long sequence of
events.

  * •

Our neural network approach can use long sequences of events with hundreds of
thousands of events.

  * •

We apply our approach to the problem of macroeconomic indexes. We predict the
daily portfolio default rate using financial transactions events of bank
clients.

  * •

Our models work better than alternatives including classic machine learning
models with hand-crafted features.

  * •

Obtained day embeddings are meaningful, and thus can be used to solve other
downstream tasks.

Figure 1: Transactions processing scheme for our approach: from an ocean of
available client financial transactions made during a particular day we sample
a useful part, then we process them via a neural network to produce a
prediction of a macroeconomic indexes for these day. The scheme allows end2end
learning of neural network parameters to directly predict the target, which is
the portfolio default rate in our case.

##  II Literature review

Our problem has two essential components: the problem of macroeconomic indexes
prediction without specification of particular input data and the problem of
representation learning for structure data, in particular usage of long
sequences of events as input to a deep learning model. In the review we
separate these two directions.

Economics papers consider prediction of macroeconomic indexes using external
data, in particular transactions aggregates. In [6] authors focus on now-
casting and forecasting quarterly private consumption. The authors used
several external data sources, in particular the data collected from automated
teller machines (ATMs), encompassing cash withdrawals at ATM terminals, and
points-of-sale (POS) payments with debit and credit cards. The received data
is transformed into aggregates of transactions (amount feature). In addition,
the considered data included the internet search patterns provided by Google
Trends to construct indicators of consumption behavior, traditional monthly
indicators as employment growth, index of retail sales and other. As a
benchmark authors used random-walk model and compare it with autoregressive
models (AR) built on different indicators. Detailed analysis demonstrated that
transaction indicators and traditional indexes seem to dominate the rest
indicators. Moreover, the model that used transaction aggregates and
traditional indexes as inputs is better, than the model that used only indexes
in terms of RMSE for different prediction horizon. In [7] the authors used the
debit card data to predict quarterly Norwegian total household consumption.
Aggregated debit card data come from the Norwegian retail clearing institution
and are sum of the values of all transactions within a week. As a benchmark
authors used autoregressive model built on past values of the target variable
and compare it with AR models built on other quarter predictors and with
mixed-data sampling (MIDAS) regressions built on available weekly data [8].
MIDAS regressions with debit card transactions data improve predicted accuracy
compare to standard AR models that use alternative quarter predictors over the
sample 2011Q4-2020Q1. In [9] the authors construct a graph of financial
transactions, where the edges are described with the following data:
transaction amount and frequency, which characterize the committed
transaction. The nodes in these graphs contain information about account
balance. The authors focused on credit card transactions, which form a
bipartite graph between holders of credit cards and the merchants where they
shop at. The main challenge is that these graphs are highly dynamic. Authors
applied Graph Neural Networks to encode high-dimensional data in low-
dimensional latent space. Calculating shifts of embeddings in the latent space
lead to quantitative estimation of impact of economic events (in particular
COVID-19) happening between these two dates where shifts are measured.

Representation learning is a major research area in machine learning and deep
learning and there are a lot of articles devoted to it. The goal is to
construct small-dimensional representations of complex structured data [10],
e.g. images [11] or sequences of events [12]. We will start from simple
approach moving toward more complicated ones.

In [13] the authors suggest sentence embedding method, improving simple
averaging of words in a sentence. The algorithm has two steps: (1) calculate
the weighted average of the words’ embeddings in the sentence, (2) remove the
projections of the average vectors on their first singular vector. The weight
of each word $w$ is $a/(a+p(w))$, with $a$ being a hyperparameter of the
algorithm and $p(w)$ is the word frequency, which can be estimated in several
ways. For text similarity problems this method shows better results than the
simple average of the vectors of words in a sentence, and even surpasses RNN
and LSTM models in most problems. This method is easy to implement,
unsupervised and can work with no limits in length of sequence lengths. Also
one can apply this method to the calculation embeddings from transaction data.
Here a sequence of transactions is a sentence, and a single transaction is a
word from a finite dictionary. The word frequency $p(w)$ is the frequency of
specific Merchant Category Code (MCC-code) for transaction in a dictionary.

In [14] the authors develops a model for the construction of sentence
embedding using Recurrent Neural Network (RNN) [15] with Long Short-Term
Memory (LSTM) cells. The proposed model creates semantic vector from sentence
by sequentially processing each word, extracting information that contains in
it and union all information from all words. The model accumulates a large
amount of information to the last word in sentence due to LSTM sells,
integrated in model, which allow to capture long term memory. To the last word
hidden layer of this network contains information about semantic meaning of
the sentence. RNNs have been used to process sequential data and compute their
embeddings. It worked well for quite short sentences. But when the sequence
length increases, the difficulty in training this network arises due to
gradient vanishing and exploding problems [16]. For RNNs learning long-term
patterns is quite difficult. To tackle this problems, new architectures were
developed: Long Short-Term Memory (LSTM) [17] and Gated Recurrent Unit (GRU)
[18]. However, the length of sequences on which these networks could be
trained is still limited by about 1000 tokens or ”words” due to gradient decay
over layers. That’s why constructing a good trainable neural network for long
sequences is a challenging task.

In [19] the authors claim that the problems with RNN and LSTM described above
can be dealed with, if we change the architecture of RNN a bit. The proposed
new architecture is called Independent Recurrent Neural Network (IndRNN). The
difference between this architecture from standard RNN is in independence of
neurons in the same layer, but it doesn’t mean that there is no connection
between neurons at all: neurons are connected only across layers.

The authors have shown that the IndRNN network can easily be trained on long
sequences, it is not exposed to gradient exploding and vanishing problems. In
addition, this network can be trained robustly with ReLU (Rectified Linear
Unit) activation function. To construct deeper neural network than existing
RNNs, authors propose to stack Multiple IndRNNs into one big model.
Experiments showed that proposed model with 21 layers can be trained robustly
on sequences of length over 5000. Comparison of IndRNNs with RNN and LSTM on
various problems showed superiority of proposed IndRNNs. Better performances
have been achieved on various tasks by using IndRNNs compared with the
traditional RNN and LSTM.

In [20] the authors suggest to use Convolution Neural Networks (CNN) to
calculate embeddings via solving time series classification problem, which is
quite challenging due to the specific nature of time series data: high feature
dimensionality, big data size and noisy nature of data itself. Traditionally
feature-based approaches were used to tackle this problem, but it’s obvious
that generating hand-crafted features from time-series leads to reducing the
original information. The deep learning techniques are explored to improve the
performance of traditional feature-based approaches. Convolution neural
network framework is proposed for this classification problem in order to
improve the performance of traditional methods. The difference from other
traditional approaches is in ability of CNN to automatically extract features
from time series data by using convolution and pooling operations which fully
describes the structure of the data. The experimental results for variety of
tasks in time series classification show that the proposed model outperforms
state-of-the-art methods in terms of the classification accuracy and noise
tolerance. The main advantage of CNNs compared to RNN and their derivatives is
that it can be applied to very long sequences (sequences of length that will
fit GPU memory and it’s about hundreds of thousands). The problem is in such
long sequences sequential nature of data will be lost yielding more pattern-
like nature to step forward.

One more approach that should be mentioned is using transformers for embedding
learning [21]. In this work authors present method for unsupervised embedding
learning from unlabeled data. A new state-of-the-art unsupervised method based
on pre-trained Transformers and Sequential Denoising Auto-Encoder (TSDAE)
which outperforms previous approaches by up to 6.4 points. It can achieve up
to 93.1% of the performance of in-domain supervised approaches. Further,
authors show that TSDAE is a strong pre-training method for learning sentence
embeddings, significantly outperforming other approaches like Masked Language
Model. A crucial shortcoming of previous studies is the narrow evaluation:
most work mainly evaluates on the single task of Semantic Textual Similarity
(STS), which does not require any domain knowledge. It is unclear if these
proposed methods generalize to other domains and tasks. Authors fill this gap
and evaluate TSDAE and other recent approaches on four different datasets from
heterogeneous domains. The main disadvantage of using transformers in their
inability to process long sequences due to their self-attention operation,
which scales quadratically with the sequence length. So, straightforward
application of neural networks to long sequences is challenging with even more
difficulties for architectures like transformers [22].

We can conclude, that the supervised prediction of macro indexes by an end2end
approach based on raw financial transaction data is absent in the literature.
The main challenge here is how should we construct embeddings for such long
sequences. We propose to conduct a research to answer these questions and
evaluate feasibility of such ideas.

##  III Methods

###  III-A Data structure

We consider two datasets: the one available to the bank on prediction of daily
credit default rate for a portfolio of clients and the open one related to the
prediction of a weekday for financial transactions. The overall
characteristics of datasets are in Table I.

During numerical experiments we split the data, so the training sample
includes 80% of the data, and the remaining 20% is for the validation. The
validation data are later in time, than the training data.

Feature | Day of Week | Default rate
---|---|---
Target | 7 classes | Float
MCC-code | 184 | 396
Transaction type | 77 | 134
Currency | - | 111
Country | - | 303
Time | Number of hours from 00:00
Amount | Float
TABLE I: Features of datasets marked by their target variable. For categorical
variable we provide number of unique values presented in the corresponding
dataset.

#### Day-of-Week data.

To provide wider understanding of strengths and weaknesses of our methods we
as well consider open transactions data on prediction of a day-of-week on the
base of clients’ transactions.

The open dataset [23] consists of the history of transactions of bank clients
for one year and three months. Information about day of week for these
transactions is available for each day. Number of transactions in dataset is 7
millions, with maximum number of transactions per day 20500, minimal number is
5600, mean number is 15000. The data include 4 input features for each
transaction and target labels presented in the above table.

Our idea was to check if we can predict day of week by transactions committed
during this day.

#### Default rate data.

The source of data for current research is a large Russian bank. The data
consists of more than $1$ billion transactions records for $2$ years and
overall portfolio default rate for each day from this period. From this
dataset $110$ millions of transactions were uniformly randomly selected for
the research to speed up experiments, as further dataset increase doesn’t lead
to the model improvement. For each day we have about 130 thousands of
transactions, making straightforward application of deep learning methods
impossible. The data includes 6 input features for each transaction and the
target portfolio default rate labels.

###  III-B Predictive model

There are tremendous amount of available data, so we should deal with them in
an effective manner.

The idea is to use a standard neural network architecture for the processing
of long transactions sequences and apply smart sampling technique during
training, which is described in section below. This method will allow to
operate all amount of available data. Neural network here play a role of
feature extractor: by ejecting useful patterns from long sequential data, we
will construct embedding, which can be used either in classification task or
regression task.

Specific details of experiments, including neural network parameters are
providen in the corresponding section.

###  III-C Training

The proposed training scheme is in Figure 2. As we can’t use all transactions
during an epoch, as they don’t fit into available memory, we sample a sequence
of uniformly random $N$ different sequential transactions. Thus, during the
whole training we’ll use the whole dataset. The training algorithm has the
following steps:

  1. 1.

Sample a sequence of uniformly random $N$ different sequential transactions
for each day from the whole dataset.

  2. 2.

Use embedding layers with learnable weights for categorical variables. These
layers map every token of its type (MCC-code/transaction
type/country/currency) to a vector space, with embedding dimension as a
parameter.

  3. 3.

Concatenate $N$ tokens into one sequence.

  4. 4.

The resulting sequences serves as the input to a neural network
(LSTM/CNN/IndRNN), which output is followed with a pooling layer and a fully
connected layer.

Figure 2: Proposed scheme for the training of a neural network model

###  III-D Inference

During the inference procedure for the classification task we sample
transactions several times from dataset, feed it to neural network, average
embeddings and make predictions on obtained embeddings.

During the inference procedure for the regression task we make predictions on
one subset of transactions and then apply moving average on predicted values.

##  IV Experiments

###  IV-A Methods details

#### Baseline.

For every calendar day we created handcrafted statistical features from
transaction data. We calculated frequencies of each MCC-code, type, country
and currency types per day and features, characterizing amount distribution.
These features were obtained by grouping amount distribution into 10 bins and
calculating mean value of every bin and and number of observations in it. So
each day was represented with the 965 tabular features described above. Via
XGBoost library [24] we train a gradient boosting with these features as
input. The best hyperparameters found via grid search for XGBoost are the
following: number of trees in ensemble is 50, depth of tree is 3, minimum sum
of instance weight (Hessian) needed in a child is 1.

#### Word2vec baseline.

We used raw transaction data with preprocessed categorical variables (MCC-
code, transaction type) via Word2vec [25]. A good baseline sequence embedding
method [13] provided unsupervised high-dimensional embeddings. On top of these
embeddings we use gradient boosting algorithm to predict the target [24].

In our case, word frequency is the frequency of MCC-code. The value of
hyperparameter in the method was selected according to best metric obtained
with boosting model on computed embeddings. Influence of $a$ parameter on
embedding space was additionally investigated. The best hyperparameters found
via grid search for XGBoost are the following: number of trees in ensemble is
75, depth of tree is 4, minimum sum of instance weight (Hessian) needed in a
child is 1.

#### CNN.

We used raw transaction data with preprocessed categorical variables (MCC-
code, transaction type) via embedding layer with learnable weights.
Experiments were conducted with sequences of lengths from 500 to 5000 with a
step of 500. Every epoch we resample sequence with the length of the specified
above from dataset. Train includes 80% of data, and 20% is for the validation.
During validation procedure we sample sequences for 30 times from dataset,
return embeddings and evaluate model accuracy on the average of the
embeddings. The number of training epochs is $500$, initial learning rate is
$\alpha=0.01$ with step decay coefficient of 0.95 and step size of 5. The
following number of training epochs was chosen to provide the best result on
validation set. We use Adam optimizer [26]. Weights are initialized randomly.
Embedding sizes are following: $77$, $25$ for MCC-code and transaction type
respectively. The batch size is $10$.

#### IndRNN.

The difference between IndRNN from standard RNN in independence of neurons in
the same layer, but it doesn’t mean that there is no connection between
neurons at all: neurons are connected only across layers.
Experiments were conducted with sequences of the following lengths: 3000,
5000, 10000. Every epoch we resample sequence with a specified length from the
dataset. The number of training epochs is $1000$, initial learning rate is
$\alpha=0.001$ with step decay coefficient of 0.95 and step size of 5. The
following number of training epochs was chosen to provide the best result on
validation set. We use standard Adam optimizer. Weights are initialized
randomly. Embedding sizes are following $36$, $18$, $16$, $16$ for MCC-code,
transaction type, country and currency respectively. Hidden size is $100$.
Batch size is $4$.

#### LSTM.

We used raw transaction data with preprocessed categorical variables MCC-code
and transaction type via an embedding layer with learnable embeddings matrix.
We consider two different types of loss-functions: a first standard cross-
entropy loss and a second triplet loss [27] for directly optimizing
embeddings. Experiments were conducted with sequences of lengths from 500 to
1000. Every epoch we resampled sequence with the length of the specified above
from dataset. During validation procedure we sample sequences for 30 times
from dataset, return embeddings and evaluate model accuracy on the average of
the embeddings. The number of training epochs is $100$, initial learning rate
is $\alpha=0.0001$. The following number of training epochs was chosen to
provide the best result on validation set. We use standard Adam optimizer.
Weights are initialized randomly. Embedding sizes are following: $77$, $25$
for MCC-code and type respectively. Hidden size is $100$. Batch size is $15$.

#### t-SNE

One of our main goals is a construction of useful data representations. To
examine our progress toward this goal, we embed our representations via
T-distributed Stochastic Neighbor Embedding t-SNE [28]. We plot 2-dimensional
t-SNE representation of embeddings obtained via neural network models and
manually-designed features prepared as input for gradient boosting model.

###  IV-B Classification of the day of week

The first problem we consider for demonstration purposes is the classification
of a day of week based on embeddings for this day. The target accuracy is the
percentage of correctly detected days of week.

Best achieved results with its corresponded parameters are in Table II. Neural
network models work better, than baselines, making them a useful tool for
processing of huge amount of data due to representation learning ability. LSTM
neural network performance is worse: its quality is similar to Gradient
boosting with Word2vec input features. We suppose, that the reason is the high
volatility of data, as we pick random subset of transactions for a day. CNN on
transaction sequences with length 2000 provides the best accuracy among all
considered methods.

Model | Number of transactions | Accuracy
---|---|---
| sampled every epoch |
XGBoost | - | $0.45\pm 0.01$
XGboost, Word2vec features | - | $0.55\pm 0.02$
LSTM, triplet loss | 500 | $0.54\pm 0.02$
LSTM, cross-entropy loss | 500 | $0.58\pm 0.02$
CNN, cross-entropy loss | 2000 | $\mathbf{0.70}\pm\mathbf{0.02}$
TABLE II: Performance comparison for day-of-week prediction: CNN model shows
the best results

CNN provides the best accuracy, so we construct embeddings for it and for
manually-selected features. Comparison of the results obtained with simple
baseline on statistical features, NLP and CNN models on raw data are in Figure
3. As expected from accuracy metric, CNN provides much better separation.

Figure 3: t-SNE embeddings computed with three models for days from model
embeddings. Each color corresponds to a specific day-of-week. _The left_
figure presents result for the baseline, _the middle one_ presents result for
the word2vec features, _the right one_ shows obtained embeddings for the CNN
model with the best separation in the embeddings space.

###  IV-C Regression of the default rate

The two main classes of approaches for supervised regression problem for
transaction sequence data are boosting on handcrafted features and neural-
network-based models, such as LSTM and CNN. These approaches were implemented
and compared with each other.

The present results are in Table III. To speak one at a time, the gradient
boosting model performs worse than other models with negative validation
$R^{2}$ score. It shows, how much more sophisticated and complex problem is
the default rate prediction compare to day of week prediction. IndRNN model
replacing its predecessor LSTM shows better results compare to the baseline
with low positive $R^{2}$ score for test set. The amount of available GPU
memory limits the maximum sequences length for IndRNN model. CNN model show
adequate performance for default rate prediction problem. Experiments show,
that performance of best model do not depend from sequence length. Moreover,
the CNN model trains 24 faster, than the IndRNN model.

Model | N | $R^{2}$
---|---|---
Baseline XGBoost | - | $-0.82\pm{0.06}$
IndRNN | 5000 | $0.19\pm{0.01}$
CNN | 10000 | $\mathbf{0.40}\pm\mathbf{0.02}$
TABLE III: Performance comparison for the default rate prediction: CNN model
shows the best results. We also show optimal processed sequence length $N$.

###  IV-D Stability of the embeddings

Our approach as well produces embeddings for days in our data. As we claim in
previous sections, such embeddings are useful, when predicting the portfolio
default rate. We hope, that these embeddings can be useful in other problems,
and are stable enough, so one can easily retrain them using new data and
obtain similar results. The stability question arises in this work, because
data is quite noisy and one should define how obtained embeddings are
appropriate at all. What if they are random? Thats’s why we define embeddings
stability metrics, and then measure it for our embeddings.

#### Embeddings stability quality metric

Given a set $S$ of $N$ elements $S=\\{s_{1},s_{2},\ldots s_{N}\\}$, consider
two partitions of $S$, namely $U=\\{U_{1},U_{2},\ldots,U_{R}\\}$ with $R$
clusters, and $V=\\{V_{1},V_{2},\ldots,V_{C}\\}$ with $C$ clusters. Our goal
is to measure correspondence between two partitions: how similar is $U$ to
$V$. We do it via Adjusted Mutual Information (AMI) score [29].

The entropy of the partition $U$ is $H(U)$, the entropy of the partition $V$
is $H(V)$. $MI(U,V)$ is the mutual information between two partitions. Then,
Adjusted Mutual Information (AMI) is calculated in the following way:

$AMI(U,V)=\frac{MI(U,V)-\mathbb{E}\\{MI(U,V)\\}}{\max\\{H(U),H(V)\\}-\mathbb{E}\\{MI(U,V)\\}},$

where $\mathbb{E}$ corresponds to expected values for a random partition. AMI
is an adjusted partition quality metric, such that for a random partitions $U$
and $V$ it will have value $0$, and the ideal correspondence will have AMI
value $1$ with bigger values are better.

#### AMI estimation

To measure stability of embeddings we estimate AMI for our embeddings in the
following way:

  * •

Sample transactions 10 times from original dataset.

  * •

Obtain 10 sets of embeddings with a separate training of a suggested
algorithm.

  * •

Cluster them via K-means clustering with 7 clusters.

  * •

Calculate AMI for obtained between each pair from 45 possible pairs.

  * •

Average 45 obtained values.

Mean AMI is $0.33$ versus $0.00$ for random embeddings. Thus, we conclude that
our embedding vectors are stable enough to be used in other problem and with
new data.

##  V Conclusion

The current work investigates the problem of macroeconomic indexes prediction
based on financial transactions data. The challenge here is the tremendous
amount of available semi-structured transactions: millions of clients commit
transactions on a specific day.

We propose an approach based on a sampling of transaction sequences and
processing these sequences via a neural network model. The evidence of the
superiority of our approach comes from results for two problems: a portfolio
of credits default rate prediction and a day-of-week prediction.

We consider two main directions for constructing such models: classic feature-
based machine learning and representation learning via neural networks. For
the classic gradient boosting model, we use hand-crafted features and word2vec
features. For neural-network-based models LSTM and CNN, we use sampling and
end2end learning for sampled sequences of transactions. CNN provides the best
accuracy of 0.70 for day-of-week prediction. According to the results for this
problem, we picked up two models: boosting on hand-crafted features and
end2end CNN on raw data to continue with them for the default rate prediction
task. For the second problem, the best $R^{2}=0.4$ is again for CNN. So, it
seems reasonable to apply CNN for such long sequences with hundreds of
thousands of transactions.

Our neural network predicts the default rate with high quality — using only
raw transactions data. The reasons are a new sampling scheme and an end2end
training procedure. The demonstrated approach shows the correlation between
the people’s transaction activity and the portfolio default rate and predicts
the future default rate. Bank employees can use obtained embeddings to solve
other downstream tasks making them universal indicators of the current
macroeconomic state.

## Acknowledgments

The work is supported by the Russian Science Foundation (project 21-11-00373).

## References

  * [1] P. Deshpande, “How COVID-19 Has Infected AI Models,” https://www.dominodatalab.com/blog/how-covid-19-has-infected-ai-models/, 2020, [Online; accessed 05-July-2021].
  * [2] M. R. Sousa, J. Gama, and M. J. S. Gonçalves, “A two-stage model for dealing with temporal degradation of credit scoring,” _arXiv preprint arXiv:1406.7775_ , 2014.
  * [3] M. R. Sousa, J. Gama, and E. Brandão, “A new dynamic modeling framework for credit risk assessment,” _Expert Systems with Applications_ , vol. 45, pp. 341–351, 2016.
  * [4] D. Babaev, M. Savchenko, A. Tuzhilin, and D. Umerenkov, “ET-RNN: Applying deep learning to credit loan applications,” in _Proceedings of the 25th ACM SIGKDD international conference on knowledge discovery & data mining_, 2019, pp. 2183–2190.
  * [5] M. Cotugno, V. Stefanelli, and G. Torluccio, “Relationship lending, default rate and loan portfolio quality,” _Applied financial economics_ , vol. 23, no. 7, pp. 573–587, 2013.
  * [6] M. Gil, J. J. Pérez, A. J. Sanchez Fuentes, and A. Urtasun, “Nowcasting private consumption: traditional indicators, uncertainty measures, credit cards and some internet data,” 2018.
  * [7] K. A. Aastveit, T. M. Fastbø, E. Granziera, K. S. Paulsen, and K. N. Torstensen, “Nowcasting norwegian household consumption with debit card transaction data,” Norges bank research, Tech. Rep., 2020.
  * [8] E. Ghysels, S.-C. Pedro, and R. Valkanov, “The midas touch: Mixed data sampling regression models,” _UCLA: Finance_ , 2004.
  * [9] A. Gogoglou, B. Nguyen, A. Salimov, J. Rider, and C. B. Bruss, “Navigating the dynamics of financial embeddings over time,” _arXiv_ , 2020.
  * [10] I. Goodfellow, Y. Bengio, and A. Courville, _Deep learning_. MIT press, 2016.
  * [11] K. Weiss, T. M. Khoshgoftaar, and D. Wang, “A survey of transfer learning,” _Journal of Big data_ , vol. 3, no. 1, pp. 1–40, 2016.
  * [12] V. Zhuzhel, R. Rivera-Castro, N. Kaploukhaya, L. Mironova, A. Zaytsev, and E. Burnaev, “Cohortney: Deep clustering for heterogeneous event sequences,” _arXiv preprint arXiv:2104.01440_ , 2021.
  * [13] S. Arora, Y. Liang, and T. Ma, “A simple but tough-to-beat baseline for sentence embeddings,” _ICLR_ , 2017.
  * [14] H. Palangi, L. Deng, Y. Shen, J. Gao, X. He, J. Chen, X. Song, and R. Ward, “Deep sentence embedding using long short-term memory networks: Analysis and application to information retrieval,” _arXiv_ , 2016.
  * [15] A. Sherstinsky, “Fundamentals of recurrent neural network (rnn) and long short-term memory (lstm) network,” _Physica D: Nonlinear Phenomena_ , vol. 404, p. 132306, Mar 2020. [Online]. Available: http://dx.doi.org/10.1016/j.physd.2019.132306
  * [16] A. Kag, Z. Zhang, and V. Saligrama, “RNNs incrementally evolving on an equilibrium manifold: A panacea for vanishing and exploding gradients?” in _International Conference on Learning Representations_ , 2019.
  * [17] S. Hochreiter and J. Schmidhuber, “Long short-term memory,” _Neural computation_ , vol. 9, no. 8, pp. 1735–1780, 1997.
  * [18] K. Cho, B. Van Merriënboer, C. Gulcehre, D. Bahdanau, F. Bougares, H. Schwenk, and Y. Bengio, “Learning phrase representations using rnn encoder-decoder for statistical machine translation,” _arXiv preprint arXiv:1406.1078_ , 2014.
  * [19] S. Li, W. Li, C. Cook, C. Zhu, and Y. Gao, “Independently recurrent neural network (indrnn): Building a longer and deeper rnn,” _arXiv_ , 2018.
  * [20] B. Zhao, H. Lu, S. Chen, J. Liu, and D. Wu, “Convolutional neural networks for time series classification,” _Journal of Systems Engineering and Electronics_ , 2017.
  * [21] K. Wang, N. Reimers, and I. Gurevych, “Tsdae: Using transformer-based sequential denoising auto-encoder for unsupervised sentence embedding learning,” _arXiv_ , 2021.
  * [22] Y. Tay, M. Dehghani, S. Abnar, Y. Shen, D. Bahri, P. Pham, J. Rao, L. Yang, S. Ruder, and D. Metzler, “Long range arena: A benchmark for efficient transformers,” _arXiv preprint arXiv:2011.04006_ , 2020.
  * [23] I. Fursov, M. Morozov, N. Kaploukhaya, E. Kovtun, R. Rivera-Castro, G. Gusev, D. Babaev, I. Kireev, A. Zaytsev, and E. Burnaev, “Adversarial attacks on deep models for financial transaction records,” _arXiv preprint arXiv:2106.08361_ , 2021.
  * [24] T. Chen and C. Guestrin, “Xgboost: A scalable tree boosting system,” in _Proceedings of the 22nd acm sigkdd international conference on knowledge discovery and data mining_ , 2016, pp. 785–794.
  * [25] T. Mikolov, K. Chen, G. Corrado, and J. Dean, “Efficient estimation of word representations in vector space,” _arXiv preprint arXiv:1301.3781_ , 2013\.
  * [26] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” _arXiv preprint arXiv:1412.6980_ , 2014.
  * [27] E. Hoffer and N. Ailon, “Deep metric learning using triplet network,” _arXiv_ , 2014.
  * [28] L. Van der Maaten and G. Hinton, “Visualizing data using t-SNE.” _Journal of machine learning research_ , vol. 9, no. 11, 2008.
  * [29] N. X. Vinh, J. Epps, and J. Bailey, “Information theoretic measures for clusterings comparison: Variants, properties, normalization and correction for chance,” _The Journal of Machine Learning Research_ , vol. 11, pp. 2837–2854, 2010.


Jinzhu Jia<EMAIL_ADDRESS>
# Non-linear Mendelian randomization with Two-stage prediction estimation and
Control function estimation

Xinpei Wang11affiliationmark:  Tao Huang2,3,42,3,4affiliationmark:  and Jinzhu
Jia1,5*1,5*affiliationmark:  11affiliationmark:  Department of Biostatistics,
School of Public Health, Peking University, 100191 China.
22affiliationmark:  Department of Epidemiology & Biostatistics, School of
Public Health, Peking University, 100191 China.
33affiliationmark:  Center for Intelligent Public Health, Academy for
Artificial Intelligence, Peking University, Beijing, China.
44affiliationmark:  Key Laboratory of Molecular Cardiovascular Sciences
(Peking University), Ministry of Education, 100191 China.
55affiliationmark:  Center for Statistical Science, Peking University, 100871
China.

###### Abstract

Most of the existing Mendelian randomization (MR) methods are limited by the
assumption of linear causality between exposure and outcome, and the
development of new non-linear MR methods is highly desirable. We introduce
two-stage prediction estimation and control function estimation from
econometrics to MR and extend them to non-linear causality. We give conditions
for parameter identification and theoretically prove the consistency and
asymptotic normality of the estimates. We compare the two methods
theoretically under both linear and non-linear causality. We also extend the
control function estimation to a more flexible semi-parametric framework
without detailed parametric specifications of causality. Extensive simulations
numerically corroborate our theoretical results. Application to UK Biobank
data reveals non-linear causal relationships between sleep duration and
systolic/diastolic blood pressure.

###### keywords:

non-linear Mendelian randomization; instrumental variable

##  1 Introduction

The presence of unmeasured confounders in epidemiological studies usually
presents a major bias for the estimation of causal effects between exposure
and outcome. Mendelian randomization (MR) uses genetic variation as an
instrumental variable (IV) to estimate the causal effect of exposure on
outcome through the association of genetic variation with exposure and
outcome, respectively Sekula et al. (2016). It is effective in avoiding
confounding bias and is now widely used in epidemiological studies. The most
commonly used MR methods are the ratio method as well as the two-stage least
squares method, both of which rely on the linear assumption of causal effects,
i.e. that each unit change in exposure produces the same effect on the outcome
Burgess et al. (2017). In practice, however, this assumption may not be
satisfied in many cases. For example, cohort study and meta-analysis have
found obvious U-shaped correlations between sleep duration and several
cardiovascular events Yin et al. (2017); Daghlas et al. (2019). Therefore,
non-linear MR methods are essential.

Researchers have already noted this problem. For example, Burgess Burgess et
al. (2014); Staley and Burgess (2017); Coscia et al. (2022); Mason and Burgess
(2022); Tian et al. (2023), as well as Silverwood Silverwood et al. (2014),
have suggested that it is possible to stratify the individuals in a certain
way, treat the causal effects as linear within each stratum, use the ratio
method to calculate local average causal effect within each stratum, and then
portray the non-linear causal effect of exposure on the outcome by testing for
heterogeneity among the stratum estimates, conducting trend tests, and fitting
fractional polynomial model or piecewise linear model. Using these methods,
epidemiological researchers have found many non-linear causal effects of
exposure on outcomes Sun et al. (2019); Rogne et al. (2020); Ai et al. (2021).
Sulc et al. Sulc et al. (2022) instead used polynomials to model the non-
linear causal relationship between exposure and outcome, again using a two-
stage approach to estimate the coefficients of each polynomial term. While
these methods have been widely used in applications, the theoretical
properties of the estimators have not been verified.

The presence of unmeasured confounders is likewise an important issue and is
referred to as endogeneity in economic studies. Economists have also developed
a number of methods to deal with unmeasured confounders. The most commonly
used of these are two-stage prediction estimation and control function
estimation Horowitz (2011); Terza et al. (2008). Both methods exhibit good
performance under linear causality but differs when extended to non-linear
causality Guo and Small (2016); Terza et al. (2008).

In this article, we introduce two-stage prediction estimation and control
function estimation from economics to MR and further refine and extend them.
We give conditions for parameter identification under non-linear causality for
both methods, theoretically prove the consistency and asymptotic normality of
the estimators, compare the two estimates under linear and non-linear
causality, and extend the control function estimation to a more flexible semi-
parametric framework. We validated our results numerically using a large
number of simulation experiments and applied the methods to UK Biobank (UKB)
data to investigate the causal relationship between sleep duration and
systolic/diastolic blood pressure.

The paper is organized as follows: in Section 2 we give the parameter
identification conditions for the two methods, in Section 3 we prove the
consistency and asymptotic normality of the estimates and compare the two
methods in both linear and non-linear causality, in Section 4 we give a more
flexible framework for semi-parametric estimation, in Section 5 we perform
numerical simulation experiments, in Section 6 we apply the methods to UKB
data, and Section 7 and Section 8 contain our discussion and conclusions,
respectively. All proofs are placed in the Appendix.

##  2 Identification

###  2.1 Setup and notation

We have i.i.d. observations $\\{(Z_{i},C_{i},X_{i},Y_{i})\\}_{i=1}^{n}$ from a
super-population $(Z,C,X,Y)$ where $Z$ is the IV, $C$ is the measured
covariate, $X$ is the exposure variable, $Y$ is the outcome, $U$ is the
unmeasured confounder. For IV estimation, $Z$ should meet three core
assumptions: (1) $Z$ has an effect on the exposure $X$; (2) $Z$ is independent
of the unmeasured confounder $U$; (3) $Z$ is independent of outcome $Y$ given
$X,C$ and $U$. In addition, we also assume $Z\Vbar C$ and $C\Vbar U$. The
relationship could be depicted using Figure 1.

Figure 1: The Causal Graph for $(Z,C,U,X,Y)$ with $U$ the unmeasured
confounder.

Considering the the most prevalent additive genetic model, it is natural to
assume that the effect of $Z$ on $X$ is linear. Without loss of generality, we
also assume that the causal effects of $C,U$ on $X$ is linear, but the causal
effects on $Y$ can be non-linear. Distinguishing from linear MR, we consider
the case where $X$ has a non-linear causal effect on $Y$, denoting as $f(X)$.
Thus, we have the following equations:

$X=\beta_{X0}+\beta_{Z}Z+\beta_{C}C+\beta_{U}U+\epsilon_{X}$ (1)
$Y=\beta_{Y0}+f(X)+g(C)+h(U)+\epsilon_{Y}$ (2)

where $\epsilon_{X}$ and $\epsilon_{Y}$ are the random error terms with
$E(\epsilon_{X})=0,E(\epsilon_{Y})=0$. $f(X)$ and $g(C)$ can be composed of
several different basis functions, i.e.,
$f(X)=\sum_{j=1}^{K_{1}}\theta_{j}f_{j}(X)$,
$g(C)=\sum_{j=1}^{K_{2}}\gamma_{j}g_{j}(C)$. Matrix
$F=[1,f_{1}(X),\ldots,f_{K_{1}}(X)]$ and matrix
$G=[1,g_{1}(C),\ldots,g_{K_{2}}(C)]$ are both full column rank. There can also
be multiple IVs, but it is required that these IVs are linearly independent,
and the measeured covariates are the same.

Then we can estimate
$\beta_{Y0},\theta_{1},\ldots,\theta_{K_{1}},\gamma_{1},\ldots,\gamma_{K_{2}}$
using the following two methods.

###  2.2 Two-stage prediction estimation

Rewrite equation 2 as:

$\begin{split}Y&=\beta_{Y0}+\sum_{j=1}^{K_{1}}\theta_{j}\hat{f}_{j}(X)+\sum_{j=1}^{K_{2}}\gamma_{j}g_{j}(C)+\\\
&\quad\sum_{j=1}^{K_{1}}\theta_{j}\left(f_{j}(X)-\hat{f}_{j}(X)\right)+h(U)+\epsilon_{Y}\end{split}$
(3)

where $\hat{f}_{j}(X)$ is the fitted values of linear regression $f_{j}(X)\sim
Z+C,j=1,2,\ldots,K_{1}$. It can be proved that

$E(Y|\hat{f}(X),g(C))=\beta_{Y0}+\sum_{j=1}^{K_{1}}\theta_{j}\hat{f}_{j}(X)+\sum_{j=1}^{K_{2}}\gamma_{j}g_{j}(C)$

So we can fit the regression of
$Y\sim\hat{f}_{1}(X)+\hat{f}_{2}(X)+\ldots+\hat{f}_{K_{1}}(X)+g_{1}(C)+g_{2}(C)+\ldots+g_{K_{2}}(C)$
to get the estimates of
$\beta_{Y0},\theta_{1},\ldots,\theta_{K_{1}},\gamma_{1},\ldots,\gamma_{K_{2}}$.

######  Theorem 1.

Denote the number of IVs and the number of measured covariates as $n_{1}$ and
$n_{2}$, respectively. The number of linear functions of $C$ in
$g_{1}(C),g_{2}(C),\ldots,g_{k2}(C)$ is $k$. Then, the parameters
$\beta_{Y0},\theta_{1},\ldots,\theta_{K_{1}},\gamma_{1},\ldots,\gamma_{K_{2}}$
can be identified if and only if $K_{1}+k\leq(n_{1}+n_{2})$.

The proof of Theorem 1 is placed in Appendix A.

######  Remark 1.

When $f(X)=X$ as well as $g(C)=C$, the two-stage prediction estimation
degenerates into the ordinary two-stage least squares estimation. Different
from the parameter identification condition $(K_{1}\leq n_{1})$ of the
ordinary two-stage least squares estimation, we allow for the non-linear
causality of both $X$ and $C$ on $Y$ and take the form of causality (linear or
non-linear) between $C$ and $Y$ into account in the parameter identification
condition.

###  2.3 Control function estimation

Rewrite equation 1 and 2 as:

$X=\beta_{X0}+\beta_{Z}Z+\beta_{C}C+\delta_{1}$ (4)
$Y=\beta_{Y0}+f(X)+g(C)+\delta_{2}$ (5)

where $\delta_{1}=\beta_{U}U+\epsilon_{X},\delta_{2}=h(U)+\epsilon_{Y}$.
Assume that the relationship between $\delta_{2}$ and $\delta_{1}$ is linear,
i.e. $\delta_{2}=\rho\delta_{1}+e,\delta_{1}\Vbar e$. Then we have

$Y=\beta_{Y0}+f(X)+g(C)+\rho\hat{\delta}_{1}+\rho(\delta_{1}-\hat{\delta}_{1})+e$
(6)

where $\hat{\delta}_{1}$ is the residuals of the linear regression $X\sim
Z+C$. It can be proved that

$E(Y|f(X),g(C),\hat{\delta}_{1})=\beta_{Y0}+f(X)+g(C)+\rho\hat{\delta}_{1}$

So we can fit the regression of $Y\sim
f_{1}(X)+\ldots+f_{K_{1}}(X)+g_{1}(C)+\ldots+g_{K_{2}}(C)+\hat{\delta}_{1}$ to
get the estimates of
$\beta_{Y0},\theta_{1},\ldots,\theta_{K_{1}},\gamma_{1},\ldots,\gamma_{K_{2}}$.

######  Theorem 2.

The parameters
$\beta_{Y0},\theta_{1},\ldots,\theta_{K_{1}},\gamma_{1},\ldots,\gamma_{K_{2}}$
can be identified when $\delta_{2}=\rho\delta_{1}+e(\delta_{1}\Vbar e)$.

The proof of Theorem 2 is placed in Appendix B.

######  Remark 2.

The parameter identification condition of the control function estimation has
been stated by several researchers in economics Guo and Small (2016);
Wooldridge (2015). We state that, different from two-stage prediction
estimation, the control function estimation has no requirements on the number
of IVs, measures covariates and the form of causality (linear or non-linear)
between $C$ and $Y$.

##  3 Estimation

###  3.1 Two-stage prediction estimation

In the first stage, we fit the regression of $f_{j}(X)\sim Z+C$ and get the
fitted values $\hat{f}_{j}(X),j=1,2,\ldots,K_{1}$. Then in the second stage,
we fit the regression of
$Y\sim\hat{f}_{1}(X)+\hat{f}_{2}(X)+\ldots+\hat{f}_{K_{1}}(X)+g_{1}(C)+g_{2}(C)+\ldots+g_{K_{2}}(C)$
to get the estimates of
$\beta_{Y0},\theta_{1},\theta_{2},\ldots,\theta_{K_{1}},\gamma_{1},\gamma_{2},\ldots,\gamma_{K_{2}}$,
denoting as $\hat{B}_{2SP}$. And the true value of parameters
$\beta_{Y0},\theta_{1},\theta_{2},\ldots,\theta_{K_{1}},\gamma_{1},\gamma_{2},\ldots,\gamma_{K_{2}}$
is denoted as $B$.

######  Theorem 3.

Under the assumption that: $K_{1}+k\leq(n_{1}+n_{2})$ and $g(C)$ is a linear
function of C, we have

$\hat{B}_{2SP}\stackrel{{\scriptstyle p}}{{\rightarrow}}B$
$\sqrt{n}(\hat{B}_{2SP}-B)\stackrel{{\scriptstyle
d}}{{\rightarrow}}N(0,\Sigma_{2SP})$

where $\Sigma_{2SP}$ is the asymptotic variance of $\sqrt{n}(\hat{B}_{2SP}-B)$
and is equal to
$E(V^{\prime}V)^{-1}Var(\delta_{2}),V=[1,\hat{f}_{1},\hat{f}_{2},\ldots,\hat{f}_{K_{1}},g_{1},g_{2},\ldots,g_{K_{2}}]$.

The assumption of $K_{1}+k\leq(n_{1}+n_{2})$ is to ensure that the parameters
are identifiable. The proof of Theorem 3 is placed in Appendix C.

######  Remark 3.

Distinguishing ourselves from other studies on two-stage prediction
estimation, we consider the effect of the form of causality of $C$ on $Y$ on
obtaining a consistent estimate and give the exact formula for the asymptotic
variance in the case where the causality between $X$ and $Y$ is non-linear.
According to Theorem 3, we can obtain a consistent estimate of $B$ and the
corresponding asymptotic variance using observed data, so that we can perform
the hypothesis test $H_{0}:B=0$ to test whether the non-linear causality is
statistically significant directly.

###  3.2 Control function estimation

In the first stage, we fit the regression of $X\sim Z+C$ and get the residuals
$\hat{\delta}_{1}$. Then in the second stage, we fit the regression of $Y\sim
f_{1}(X)+f_{2}(X)+\ldots+f_{K_{1}}(X)+g_{1}(C)+g_{2}(C)+\ldots+g_{K_{2}}(C)+\hat{\delta}_{1}$
to get the estimates of
$\beta_{Y0},\theta_{1},\theta_{2},\ldots,\theta_{K_{1}},\gamma_{1},\gamma_{2},\ldots,\gamma_{K_{2}}$,
denoting as $\hat{B}_{CF}$.

######  Theorem 4.

Under the assumption that $\delta_{2}=\rho\delta_{1}+e(\delta_{1}\Vbar e)$, we
have

$\hat{B}_{CF}\stackrel{{\scriptstyle p}}{{\rightarrow}}B$
$\sqrt{n}(\hat{B}_{CF}-B)\stackrel{{\scriptstyle
d}}{{\rightarrow}}N(0,\Sigma_{CF})$

where $\Sigma_{CF}$ is the asymptotic variance of $\sqrt{n}(\hat{B}_{CF}-B)$.
(The exact form of $\Sigma_{CF}$ is somewhat complicated, and we place it in
Appendix D.2.)

The proof of Theorem 4 is placed in Appendix D.

######  Remark 4.

Similar to Theorem 3, according to Theorem 4, we can obtain a consistent
estimate of $B$ and the corresponding asymptotic variance using observed data
by control function estimation, so that we can perform the hypothesis test
$H_{0}:B=0$ to test whether the non-linear causality is statistically
significant directly.

###  3.3 Comparison of the two methods

######  Theorem 5.

On the basis that the assumptions of Theorems 1-4 are satisfied, the control
function estimation is equivalent to the two-stage prediction estimation when
the causality between $X$ and $Y$ is linear, and the control function
estimation is more efficient than the two-stage prediction estimation when the
causality between $X$ and $Y$ is non-linear.

We give the proof of Theorem 5 in the Appendix E.

######  Remark 5.

Different from Guo and Small’s Guo and Small (2016) study, we compare two-
stage prediction estimation and control function estimation from two aspects:
first, the conditions for parameter identification and consistency of
estimators (Theorem 1,2,3,4); and second, the asymptotic variance of the
estimates of the two methods when the assumptions required for both methods
are satisfied (Theorem 5).

##  4 Flexible semi-parametric MR

Considering the good performance of control function estimation when non-
linear causality is present, we can further extend it to more flexible semi-
parametric case. In concrete terms, we fit the linear regression of $X\sim
Z+C$ and get the residuals $\hat{\delta}_{1}$ just as before. But in the
second stage we adopt a non-parametric estimation, using a flexible
specification for $f(X)$, instead of the original detailed parametric
estimation.

By Theorem 4 we know that fitting the regression of $Y\sim
f_{1}(X)+f_{2}(X)+\ldots+f_{K_{1}}(X)+g(C)+\hat{\delta}_{1}$ in the second
stage yields consistent and asymptotically normal estimates of
$\theta_{j},j=1,2,\ldots,K_{1}$. Thus we can represent $f(X)$ using spline
bases and fit non-parametric penalized spline regression of $Y\sim
s(X)+s(C)+\hat{\delta}_{1}$ to estimate $f(X)$ in the second stage. Given the
property of spline smoothers Perperoglou et al. (2019), it is possible to
obtain an estimate of $f(X)$ with a very small approximation error as well as
a very small tendency to overfitting.

##  5 Simulation

We performed simulations in both the parametric and semi-parametric cases. In
the parametric case we specified the form of $f(X)$ in advance and estimated
its coefficient, while in the semi-parametric case we fitted the form of
$f(X)$ using splines. Comparisons were made between the performance of two-
stage prediction estimation and control function estimation under different
sample sizes ($N=5000,10000,20000,50000$) and different IV strengths
(proportion of exposure variance explained by IV (PVE) $=25\%,50\%,75\%$).

Figure 2: Simulation results for the parametric case. We run $1000$
simulations in each setting. The diamond represents the position of the mean
value. PVE: proportion of exposure variance explained by instrumental
variable; N: sample size; CF: control function estimation; TSP: two-stage
prediction estimation.

In both cases, $Z,C,U,\epsilon_{X},e$ were sampled from standard normal
distribution independently, $X$ and $Y$ were generated according to the
following equations, with PVE changed by varying the value of $\beta_{Z}$.

$\delta_{1}=U+\epsilon_{X}$ $\delta_{2}=\delta_{1}+e$
$X=1+\beta_{Z}Z+C+\delta_{1}$ $Y=1+f(X)+C+\delta_{2}$

Parametric Case: We set $f(X)=sin(X)$ and estimated its coefficients using
both two-stage prediction estimation and control function estimation. Figure 2
illustrates the performance of the two methods for different IV strengths and
different sample sizes (here the true true coefficient of $f(X)$ is $1$). It
can be seen that the control function method yields accurate and robust
estimates in all cases. But the two-stage prediction method shows a larger
variance, especially when the IV is stronger $(PVE=75\%)$ (Figure 2).

Semi-parametric Case: In this case, it is assumed that we do not have
information about the exact form of $f(X)$ ahead of time, the two-stage
prediction method cannot be used, and we only estimate using the control
function method. We considered four different $f(X):X,(X/3)^{2},sin(X)$, and
$e^{X/3}$. In the first stage a linear regression of $X\sim Z+C$ was fitted to
obtain the residuals $\hat{\delta}_{1}$ as in the parametric case, but in the
second stage a penalized spline regression of $Y\sim s(X)+C+\hat{\delta}_{1}$
was fitted, where $s(X)$ denotes the smooth function of $X$, the exact form of
which needs to be estimated from the data. We used cubic splines for stage
two. Figure 3 illustrates the estimates of $f(X)$ by this semi-parametric
control function method at a sample size of $5000$ and $PVE=25\%$. It is easy
to see that even with smaller sample size and weaker IV strength the method
shows good performance in estimating the specific form of $f(X)$, especially
in the denser range of the data (Figure 3).

Figure 3: Simulation results for the semi-parametric case. The title of each
subgraph is $f(X)$ set at the time of data generation.

##  6 Real data

We applied our methods to UKB real data to investigate the potential causal
effect of sleep duration on systolic as well as diastolic blood pressure. UKB
is a large biomedical database containing over 500000 UK participants Bycroft
et al. (2018). From this database, we selected White British participants,
excluding individuals identified as outliers in heterozygosity and missing
rates, excluding individuals with inconsistent genetic and self-reported sex,
excluding individuals with ten or more third-degree relatives identified and
individuals excluded from the kinship inference process, and excluding
individuals missing the exposure, outcome, and covariates we wanted to study.

We obtained 78 genetic variations significantly associated with sleep duration
from a large-scale genome-wide association study of sleep duration in a
population of European descent Dashti et al. (2019) and used unweighted
polygenic risk scores for these 78 genetic variations as IV. Mode
interpolation was used for missing genotypes. Data on sleep duration were
obtained from participants’ self-reports, and individuals with extreme data
($<3$h or $>18$h) were excluded Dashti et al. (2019). Systolic and diastolic
blood pressure were averaged using two independent automated measurements. For
individuals taking antihypertensive medication, we corrected systolic and
diastolic blood pressure by increasing 15 and 10 mmHg, respectively Evangelou
et al. (2018). In addition, we adjusted for the covariates of age, sex,
measurement center, genotyping chip, and the first 10 genetic principal
components in both stages of the regression. All of the above data were
obtained at the initial visit. A total of 402460 participants were ultimately
included in the analysis.

Figure 4: Results of real data application. SBP: systolic blood pressure; DBP:
diastolic blood pressure.

In the first stage, we fitted linear regression of sleep duration on IV as
well as covariates to obtain residuals. In the second stage, we fitted
penalized spline regressions of systolic/diastolic blood pressure on sleep
duration, residuals obtained from the first stage, and the covariates. Figure
4 gives the curves for the causal effect of sleep duration on systolic and
diastolic blood pressure derived from the use of cubic spline in the second
stage. The estimated causal effects of sleep duration on systolic and
diastolic blood pressure show a similar pattern, with changes in sleep
duration having a flat effect on systolic/diastolic blood pressure when sleep
duration is shorter than 7 hours, and prolonged sleep duration elevating
systolic/diastolic blood pressure when sleep duration is longer than 7 hours.

##  7 Discussion

Non-linear MR methods, although they have been used in many epidemiological
studies, have always lacked a theoretical proof of the mathematical properties
of the estimators. In this paper, we systematically give a theoretical study
of two-stage prediction estimation and control function estimation, including
the identification condition, consistency and asymptotic normality of the
estimates, and the difference between two methods in both linear and non-
linear causality. Overall, compared to the two-stage prediction estimation,
the control function estimation is more robust in causal parameter estimating
and more flexible in causal function form specification.

MR analysis has been widely used in epidemiologic studies in recent years due
to its good performance in dealing with confounding bias, and has gradually
developed into an important epidemiologic research method. While using MR
researchers must recognize the need for caution in interpreting its results.
Currently widely used MR methods such as inverse variance weighted method
Burgess et al. (2016), MR-Egger regression method Bowden et al. (2015),
weighted median method Bowden et al. (2016), and weighted mode method Hartwig
et al. (2017) all rely on the assumption that the causal effect between
exposure and outcome is linear. When this assumption is not satisfied, the
estimates obtained using these methods will no longer be the conditional
average treatment effect we are interested in, missing much of practical
significance. However, the causal relationship between exposure and outcome in
many cases cannot be described by a straight line simply. Thus the development
and use of MR methods applicable to non-linear causality is essential.

In terms of non-linear MR methods, a great deal of work has been done by
Burgess et al Burgess et al. (2014); Staley and Burgess (2017); Tian et al.
(2023). They have developed a series of stratification-based non-linear MR
methods and have used them to discover many non-linear causal relationships
between exposures and outcomes. In contrast to the approach of Burgess et al,
our approach does not stratify exposures, but rather tests the possible non-
linear effects of exposure on outcomes using exposure data at all levels. For
both two-stage prediction estimation and control function estimation, we give
identification and consistency conditions as well as the asymptotic variance
that allows the researchers to perform hypothesis testing on the coefficients
directly. We also compare the two methods in both linear and non-linear cases,
to facilitate the researcher’s choice of the appropriate method depending on
the specific exposure and outcome scenarios to be analyzed. In addition, we
provide a more flexible semi-parametric estimation framework, which is more
recommended if the researchers do not have any a priori knowledge about the
causal relationship to be investigated.

In the simulation experiments, under non-linear causality, we observe that the
variance of the two-stage prediction estimates is relatively large, especially
when the IV is of high strength. It is understandable that due to the linear
association between exposure and IV, the relationship between the non-linear
term of exposure and IV is inevitably non-linear, and thus the estimation
instability occurs when we use IV to linearly predict the non-linear term of
exposure. This situation is more likely to occur when the IV strength is large
or when the non-linear term of exposure deviates from linearity heavily. In
the real data application, sleep duration data were derived from participants’
self-reports and took all integer values. This leaves few values of sleep
duration less than 7 hours (only 3,4,5,6), which may be a contributing factor
to the fact that we did not observe an apparent causal effect of sleep
duration on systolic/diastolic blood pressure when sleep duration was shorter
than 7 hours.

It should be noted that our approach still has some limitations. For example,
our method is a one-sample MR method, which still needs to use individual
level data for estimation. But considering that more and more excellent large
cohorts are now opening up their data requests to researchers, we do not see
this as an insurmountable difficulty. In addition, the current two-stage
prediction estimation and control function estimation rely on a number of
assumptions, and how to obtain good estimates in relaxing these assumptions
needs for further research. Finally, in this article we have only studied the
case where both exposure and outcome are continuous variables, and estimation
when the variables are in other forms deserves to be studied further.

##  8 Conclusions

Both two-stage prediction estimation and control function estimation can yield
consistent and asymptotically normal estimates under specific conditions. The
control function estimation is much more robust than two-stage prediction
estimation under non-linear causality. In addition, the control function
estimation can be extended to a more flexible semi-parametric framework
without making detailed parametric assumptions about the form of the causal
effect.

The Authors declare that there is no conflict of interest

This work was supported by Peking University.

## References

  * Ai et al. (2021) Ai S, Zhang J, Zhao G, Wang N, Li G, So HC, Liu Y, Chau SWH, Chen J, Tan X, Jia F, Tang X, Shi J, Lu L and Wing YK (2021) Causal associations of short and long sleep durations with 12 cardiovascular diseases: linear and nonlinear mendelian randomization analyses in uk biobank.  _EUROPEAN HEART JOURNAL_ 42(34): 3349+.  10.1093/eurheartj/ehab170.
  * Bowden et al. (2015) Bowden J, Smith GD and Burgess S (2015) Mendelian randomization with invalid instruments: effect estimation and bias detection through egger regression.  _INTERNATIONAL JOURNAL OF EPIDEMIOLOGY_ 44(2): 512–525.  10.1093/ije/dyv080.
  * Bowden et al. (2016) Bowden J, Smith GD, Haycock PC and Burgess S (2016) Consistent estimation in mendelian randomization with some invalid instruments using a weighted median estimator.  _GENETIC EPIDEMIOLOGY_ 40(4): 304–314.  10.1002/gepi.21965.
  * Burgess et al. (2014) Burgess S, Davies NM, Thompson SG and Consortium EI (2014) Instrumental variable analysis with a nonlinear exposure-outcome relationship.  _EPIDEMIOLOGY_ 25(6): 877–885.  10.1097/EDE.0000000000000161.
  * Burgess et al. (2016) Burgess S, Dudbridge F and Thompson SG (2016) Combining information on multiple instrumental variables in mendelian randomization: comparison of allele score and summarized data methods.  _STATISTICS IN MEDICINE_ 35(11): 1880–1906.  10.1002/sim.6835.
  * Burgess et al. (2017) Burgess S, Small DS and Thompson SG (2017) A review of instrumental variable estimators for mendelian randomization.  _STATISTICAL METHODS IN MEDICAL RESEARCH_ 26(5, SI): 2333–2355.  10.1177/0962280215597579.
  * Bycroft et al. (2018) Bycroft C, Freeman C, Petkova D, Band G, Elliott LT, Sharp K, Motyer A, Vukcevic D, Delaneau O, O’Connell J, Cortes A, Welsh S, Young A, Effingham M, McVean G, Leslie S, Allen N, Donnelly P and Marchini J (2018) The uk biobank resource with deep phenotyping and genomic data.  _NATURE_ 562(7726): 203+.  10.1038/s41586-018-0579-z.
  * Coscia et al. (2022) Coscia C, Gill D, Benitez R, Perez T, Malats N and Burgess S (2022) Avoiding collider bias in mendelian randomization when performing stratified analyses.  _EUROPEAN JOURNAL OF EPIDEMIOLOGY_ 37(7): 671–682.  10.1007/s10654-022-00879-0.
  * Daghlas et al. (2019) Daghlas I, Dashti HS, Lane J, Aragam KG, Rutter MK, Saxena R and Vetter C (2019) Sleep duration and myocardial infarction.  _Journal of the American College of Cardiology_ 74(10): 1304—1314.  10.1016/j.jacc.2019.07.022.  URL https://europepmc.org/articles/PMC6785011.
  * Dashti et al. (2019) Dashti HS, Jones SE, Wood AR, Lane JM, van Hees VT, Wang H, Rhodes JA, Song Y, Patel K, Anderson SG, Beaumont RN, Bechtold DA, Bowden J, Cade BE, Garaulet M, Kyle SD, Little MA, Loudon AS, Luik A I, Scheer FAJL, Spiegelhalder K, Tyrrell J, Gottlieb DJ, Tiemeier H, Ray DW, Purcell SM, Frayling TM, Redline S, Lawlor DA, Rutter MK, Weedon MN and Saxena R (2019) Genome-wide association study identifies genetic loci for self-reported habitual sleep duration supported by accelerometer-derived estimates.  _NATURE COMMUNICATIONS_ 10\.  10.1038/s41467-019-08917-4.
  * Evangelou et al. (2018) Evangelou E, Warren HR, Mosen-Ansorena D, Mifsu B, Pazoki R, Gao H, Ntritsos G, Dimou N, Cabrer CP, Karaman I, Ng F, Evangelou M, Witkowska K, Tzanis E, Hellwege JN, Giri A, Edwards DRV, Sun Y V, Cho K, Gaziano JM, Wilson PWF, Tsao PS, Kovesdy CP, Esko T, Magi R, Milani L, Almgren P, Boutin T, Debette S, Ding J, Giulianini F, Holliday EG, Jackson AU, Li-Gao R, Lin WY, Luan J, Mangino M, Oldmeadow C, Prins BP, Qian Y, Sargurupremraj M, Shah N, Surendran P, Theriault S, Verweij N, Willems SM, Zhao JH, Amouyel P, Connell J, de Mutsert R, Doney ASF, Farrall M, Menni C, Morris AD, Noordam R, Pare G, Poulter NR, Shields DC, Stanton A, Thom S, Abecasis G, Amin N, Arking DE, Ayers KL, Barbieri CM, Batini C, Bis JC, Blake T, Bochud M, Boehnke M, Boerwinkle E, Boomsma D I, Bottinger EP, Braund PS, Brumat M, Campbell A, Campbell H, Chakravarti A, Chambers JC, Chauhan G, Ciullo M, Cocca M, Collins F, Cordell HJ, Davies G, de Borst MH, de Geus EJ, Deary IJ, Deelen J, Del Greco FM, Demirkale CY, Doerr M, Ehret GB, Elosua R, Enroth S, Erzurumluoglu AM, Ferreira T, Franberg M, Franco OH, Gandin I, Gasparini P, Giedraitis V, Gieger C, Girotto G, Goel A, Gow AJ, Gudnason V, Guo X, Gyllensten U, Hamsten A, Harris TB, Harris SE, Hartman CA, Havulinna AS, Hicks AA, Hofer E, Hofman A, Hottenga JJ, Huffman JE, Hwang SJ, Ingelsson E, James A, Jansen R, Jarvelin MR, Joehanes R, Johansson A, Johnson AD, Joshi PK, Jousilahti P, Jukema JW, Jula A, Kahonen M, Kathiresan S, Keavney BD, Khaw KT, Knekt P, Knight J, Kolcic I, Kooner JS, Koskinen S, Kristiansson K, Kutalik Z, Laan M, Larson M, Launer LJ, Lehne B, Lehtimaki T, Liewald DCM, Lin L, Lind L, Lindgren CM, Liu Y, Loos RJF, Lopez LM, Lu Y, Lyytikainen LP, Mahajan A, Mamasoula C, Marrugat J, Marten J, Milaneschi Y, Morgan A, Morris AP, Morrison AC, Munson PJ, Nalls MA, Nandakumar P, Nelson CP, Niiranen T, Nolte IM, Nutile T, Oldehinkel AJ, Oostra BA, O’Reilly PF, Org E, Padmanabhan S, Palmas W, Palotie A, Pattie A, Penninx BWJH, Perola M, Peters A, Polasek O, Pramstaller PP, Nguyen QT, Raitakari OT, Ren M, Rettig R, Rice K, Ridker PM, Ried JS, Riese H, Ripatti S, Robino A, Rose LM, Rotter J I, Rudan I, Ruggiero D, Saba Y, Sala CF, Salomaa V, Samani NJ, Sarin AP, Schmidt R, Schmidt H, Shrine N, Siscovick D, Smith A V, Snieder H, Sober S, Sorice R, Starr JM, Stott DJ, Strachan DP, Strawbridge RJ, Sundstrom J, Swertz MA, Taylor KD, Teumer A, Tobin MD, Tomaszewski M, Toniolo D, Traglia M, Trompet S, Tuomilehto J, Tzourio C, Uitterlinden AG, Vaez A, van der Most PJ, van Duijn CM, Vergnaud AC, Verwoert GC, Vitart V, Voelker U, Vollenweider P, Vuckovic D, Watkins H, Wild SH, Willemsen G, Wilson JF, Wright AF, Yao J, Zemunik T, Zhang W, Attia JR, Butterworth AS, Chasman D, Conen D, Cucca F, Danesh J, Hayward C, Howson JMM, Laakso M, Lakatta EG, Langenberg C, Melander O, Mook-Kanamori DO, Palmer CNA, Risch L, Scott RA, Scott RJ, Sever P, Spector TD, van der Harst P, Wareham NJ, Zeggini E, Levy D, Munroe PB, Newton-Cheh C, Brown MJ, Metspalu A, Hung AM, O’Donnell C, Edwards TL, Psaty BM, Tzoulaki I, Barnes MR, Wain L V, Elliott P, Caulfield MJ and Program MV (2018) Genetic analysis of over 1 million people identifies 535 new loci associated with blood pressure traits.  _NATURE GENETICS_ 50(10): 1412+.  10.1038/s41588-018-0205-x.
  * Guo and Small (2016) Guo Z and Small DS (2016) Control function instrumental variable estimation of nonlinear causal effect models.  _JOURNAL OF MACHINE LEARNING RESEARCH_ 17\.
  * Hartwig et al. (2017) Hartwig FP, Smith GD and Bowden J (2017) Robust inference in summary data mendelian randomization via the zero modal pleiotropy assumption.  _INTERNATIONAL JOURNAL OF EPIDEMIOLOGY_ 46(6): 1985–1998.  10.1093/ije/dyx102.
  * Horowitz (2011) Horowitz JL (2011) Applied nonparametric instrumental variables estimation.  _ECONOMETRICA_ 79(2): 347–394.  10.3982/ECTA8662.
  * Mason and Burgess (2022) Mason AM and Burgess S (2022) Software application profile: Sumnlmr, an r package that facilitates flexible and reproducible non-linear mendelian randomization analyses.  _INTERNATIONAL JOURNAL OF EPIDEMIOLOGY_ 51(6): 2014–2019.  10.1093/ije/dyac150.
  * Perperoglou et al. (2019) Perperoglou A, Sauerbrei W, Abrahamowicz M, Schmid M and Initiative TS (2019) A review of spline function procedures in r.  _BMC MEDICAL RESEARCH METHODOLOGY_ 19\.  10.1186/s12874-019-0666-3.
  * Rogne et al. (2020) Rogne T, Solligard E, Burgess S, Brumpton BM, Paulsen J, Prescott HC, Mohus RM, Gustad LT, Mehl A, Asvold BO, DeWan AT and Damas JK (2020) Body mass index and risk of dying from a bloodstream infection: A mendelian randomization study.  _PLOS MEDICINE_ 17(11).  10.1371/journal.pmed.1003413.
  * Sekula et al. (2016) Sekula P, Del Greco FM, Pattaro C and Koettgen A (2016) Mendelian randomization as an approach to assess causality using observational data.  _JOURNAL OF THE AMERICAN SOCIETY OF NEPHROLOGY_ 27(11): 3253–3265.  10.1681/ASN.2016010098.
  * Silverwood et al. (2014) Silverwood RJ, Holmes MV, Dale CE, Lawlor DA, CWhittaker J, Smith GD, Leon DA, Palmer T, Keating BJ, Zuccolo L, Casas JP, Dudbridge F and Consortium AA (2014) Testing for non-linear causal effects using a binary genotype in a mendelian randomization study: application to alcohol and cardiovascular traits.  _INTERNATIONAL JOURNAL OF EPIDEMIOLOGY_ 43(6): 1781–1790.  10.1093/ije/dyu187.
  * Staley and Burgess (2017) Staley JR and Burgess S (2017) Semiparametric methods for estimation of a nonlinear exposure-outcome relationship using instrumental variables with application to mendelian randomization.  _GENETIC EPIDEMIOLOGY_ 41(4): 341–352.  10.1002/gepi.22041.
  * Sulc et al. (2022) Sulc J, Sjaarda J and Kutalik Z (2022) Polynomial mendelian randomization reveals non-linear causal effects for obesity-related traits.  _HUMAN GENETICS AND GENOMICS ADVANCES_ 3(3).  10.1016/j.xhgg.2022.100124.
  * Sun et al. (2019) Sun YQ, Burgess S, Staley JR, Wood AM, Bell S, Kaptoge SK, Guo Q, Bolton TR, Mason AM, Butterworth AS, Di Angelantonio E, Vie GA, Bjorngaard JH, Kinge JM, Chen Y and Mai XM (2019) Body mass index and all cause mortality in hunt and uk biobank studies: linear and non-linear mendelian randomisation analyses.  _BMJ-BRITISH MEDICAL JOURNAL_ 364\.  10.1136/bmj.l1042.
  * Terza et al. (2008) Terza JV, Basu A and Rathouz PJ (2008) Two-stage residual inclusion estimation: Addressing endogeneity in health econometric modeling.  _JOURNAL OF HEALTH ECONOMICS_ 27(3): 531–543.  10.1016/j.jhealeco.2007.09.009.
  * Tian et al. (2023) Tian H, Mason AM, Liu C and Burgess S (2023) Relaxing parametric assumptions for non-linear mendelian randomization using a doubly-ranked stratification method.  _PLOS GENETICS_ 19(6).  10.1371/journal.pgen.1010823.
  * Wooldridge (2015) Wooldridge JM (2015) Control function methods in applied econometrics.  _JOURNAL OF HUMAN RESOURCES_ 50(2): 420–445.  10.3368/jhr.50.2.420.
  * Yin et al. (2017) Yin J, Jin X, Shan Z, Li S, Huang H, Li P, Peng X, Peng Z, Yu K, Bao W, Yang W, Chen X and Liu L (2017) Relationship of sleep duration with all-cause mortality and cardiovascular events: A systematic review and dose-response meta-analysis of prospective cohort studies.  _JOURNAL OF THE AMERICAN HEART ASSOCIATION_ 6(9).  10.1161/JAHA.117.005947.