footnote

# Two Stage Least Squares with Time-Varying Instruments: An Application to an
Evaluation of Treatment Intensification for Type-2 Diabetes

Daniel Tompsett  Dr Daniel Tompsett, Primary Care and Population Health,
University College London, United Kingdom, email<EMAIL_ADDRESS>Stijn
Vansteelandt  Professor Stijn Vansteelandt, Department of Applied Mathematics,
Computer Science, and Statistics, University of Ghent, Belgium, email:
<EMAIL_ADDRESS>Richard Grieve  Professor Richard Grieve, Centre
for Data and Statistical Science for Health (DASH), London School of Hygiene
and Tropical Medicine, United Kingdom, email<EMAIL_ADDRESS>Irene Petersen  Professor Irene Petersen, Primary Care and Population Health,
University College London, United Kingdom, email<EMAIL_ADDRESS>Manuel
Gomes  Professor Manuel Gomes, Primary Care and Population Health, University
College London, United Kingdom, email<EMAIL_ADDRESS>
###### Abstract

As longitudinal data becomes more available in many settings, policy makers
are increasingly interested in the effect of time-varying treatments (e.g.
sustained treatment strategies). In settings such as this, the preferred
analysis techniques are the g-methods, however these require the untestable
assumption of no unmeasured confounding. Instrumental variable analyses can
minimise bias through unmeasured confounding. Of these methods, the Two Stage
Least Squares technique is one of the most well used in Econometrics, but it
has not been fully extended, and evaluated, in full time-varying settings.
This paper proposes a robust two stage least squares method for the
econometric evaluation of time-varying treatment. Using a simulation study we
found that, unlike standard two stage least squares, it performs relatively
well across a wide range of circumstances, including model misspecification.
It compares well with recent time-varying instrument approaches via
g-estimation. We illustrate the methods in an evaluation of treatment
intensification for Type-2 Diabetes Mellitus, exploring the exogeneity in
prescribing preferences to operationalise a time-varying instrument.

Keywords: Instrumental Variable, Time-Varying, Two Stage Least Squares,
Physician Preference,Diabetes

##  1 Introduction

As routinely collected data has become more available, there has been an
increasing interest in studies identifying time-varying treatment effects in
observational databases. For example a study of diabetic treatment on blood
glucose profile over x years may be interested in both the sustained effect of
treatment as well as the trajectories of treatment effect over time. Time-
varying confounding, where, say, health related measures impact both treatment
and outcome after the first recorded time period are a major challenge. Since
time-varying confounders often serve as mediators of the effect of previous
treatment, simply controlling for them all at once blocks indirect effects of
the treatment. To handle this, sophisticated methods such as the g-methods [1]
are needed. Such methods however require that there is no unmeasured time-
varying confounding.

A popular approach to deal with unmeasured confounding is to use Instrumental
Variables (IVs), sources of exogenous variation which 1: strongly directly
affect treatment assignment, and 2: affect the outcome only through the
treatment. IV analysis has had widespread use across multiple disciplines such
as genetics, econometrics and clinical studies. Work on IV methods in single
IV, time fixed settings is numerous, with the Wald estimator, Two Stage Least
Squares (2SLS) and Inverse Variance Weighting among the most notable [2].

This has motivated research in recent years to develop methods that apply IVs
to time-varying settings. This is a uniquely difficult challenge, as to
identify multiple distinct causal effects, one requires at least as many
instruments, when identifying even one instrument can be difficult [3]. When
there are fewer IVs, strict assumptions must be applied to the causal estimand
[4]. One solution is to use a time-varying measure as an instrument.

The main challenge with applying IVs in time-varying settings is to identify
an instrument that remains strongly associated with treatment assignment over
time. Time varying IVs have recently been used in the fitting of Marginal
Structural Models using a novel inverse weighting procedure [5, 6], which used
travel distance from the nearest treatment facility as an instrument at two
recorded times. Other recent works [7] used the same time varying instrument
to improve Dynamic Treatment Regimes (DTRs) via IVs. One limitation of these
methods was relatively simple data structures with regards to the dependence
of IV on treatment and confounding. The weighting and DTR approaches only
considered a binary IV, dichotomised from a continuous measure, which tends to
reduce the strength of the IV. A recent work by the authors [8] investigated
using time varying instruments in a g-estimation setting, based on an
adaptation of recent methods [4]. It showed that a g-estimation approach may
perform better than a weighting approach using binary IVs.

The papers primary aim is to develop and investigate the application of time
varying IVs using Two Stage Least Squares (2SLS). Whilst this has been
considered [9], there is, at least to our knowledge, relatively little
guidance as to how to apply a time-varying instrument in 2SLS settings.

Of the well established IV methods, 2SLS is a strong choice. It’s easily used
and widely understood in a variety of disciplines, including genomics,
biostatistics and econometrics [10, 11, 12]. It is also a fairly modular,
comprising of two steps that can be adapted and expanded upon to meet several
challenges. Crucially, it is the basis of most work in Mendelian Randomisation
(MR) studies, of which has seen significant development of 2SLS, using genetic
traits as a rich source of multiple instruments at baseline [13, 14, 3].

We adapt existing multivariate 2SLS methods in MR [12, 15, 16] and Robust 2SLS
methods in time fixed settings [17, 18], to settings with time-varying
confounding and evaluate its statistical properties in a simulation study. We
compare the relative performance of this approach with standard 2SLS methods
and a g-estimation approach. Methods are illustrated in a case study
evaluating second line treatment for type 2 diabetes, using GP prescribing
preferences for a new class of inhibitors as the time-varying IV. Preference
is widely applied as an IV in retrospective clinical studies [11], and is
naturally time varying in nature. Whilst GP preference was considered
previously [8], the paper focused on comparison of methods, and the case study
had limited scope. We highlight the challenges of operationalising
prescription preference as a time-varying IV.

##  2 Motivating Example: Comparative Study of Second Line Therapy for Type 2
Diabetes

Our paper is motivated by an analysis of the effectiveness of second line
therapy for Type 2 Diabetes (T2D) on blood glucose levels. T2D is a
progressive disease characterized by a impaired ability for pancreatic
${\beta}$-cells to release insulin, leading to elevated glycated haemoglobin
(HbA1c), or blood glucose levels. More than 4.4 million people are estimated
to be living with T2D in the UK with more than 3.2 million at an increased
risk in the future [19], and contributes to increased risks of cardiovascular
disease, chronic kidney disease and vascular difficulties.

Treatment involves prescribed medication to control and lower HbA1c. In first
instance, NICE (2022) prescription guidelines in the UK recommend Metformin
monotherapy. However as patients grow resistant to treatment, a second line
intensification is often necessary. Second line therapy supplements Metformin
with a second oral anti-diabetic, the choice of which, while guided, is left
to clinicians and primary care practices. For this reason, second line therapy
preferences can differ greatly between practices and GPs, and is subject to
change over time [20]. Of these, patients without high risk of cardio-vascular
disease are most commonly assigned Sulfonylureas (SU), or DPP4-inhibitors
(DPP4). As one of the first intensifications for Diabetes, GPs tend to have a
strong historical preference for SU. However recent studies advocating the
advantages of DPP4 may have shifted preference towards DPP4 in recent years
[21, 22].

#### Study Population and Eligibility Criteria

The data comprises of routinely collected primary care data from three East
London clinical commissioning groups (CCGs) based in Tower Hamlets, Newham and
City, and Hackney. The data was provided by the Queen Mary University of
London, Clinical Effectiveness Group (CEG). Data on treatment and health
related information was collected and recorded in 6 monthly intervals from
2012 to 2018. Our initial data subset includes $n=7342$ patients who were
recorded as initiating second line therapy for diabetes for the first time,
after monotherapy failed, between October 2012 to October 2017. This period
coincided with a shift in the view of prescribing DPP4 versus SU [21].

The median time patients stay in the study after starting second line therapy
is 5 follow up periods, or around 2.5-3 years. We look to follow up patients
for up to 2 years, or 4 periods. Time $t=1$ is taken as the first period the
patient was recorded taking second line therapy, indicating they had begun
this treatment within the last 6 months. Times $t=2$ and $t=3$ represent a 6
and 12 months follow up time respectively, with the outcome recorded 6 months
later, 18-24 months after initiation.

Eligible patients were between 18 and 89 years of age, registered with a
primary care practice, and initiated second line therapy after first line
monotherapy failed. Patients were required to be on either SU or DPP4 at
initiation of second line therapy, with complete relevant data available for
the full follow up period. Patients who do not start on one of these two
treatment regimes, leaves the study before 3 follow up times, pauses treatment
on SU or DPP4, or begins a further intensification by taking both or another
diabetic treatment, are censored from the study. There are $n=3640$ patients
who meet this criteria with complete treatment data, and $n=2561$ with
complete data on all relevant variables.

#### Treatment Comparison

Treatment is a contrast of one of two first intensification second line
treatments.

Treatment: Initiate treatment intensification with DPP4 over 18 months.

Control: Initiate treatment intensification with SU over 18 months.

Treatment is recorded at each 6 month interval, with the ”Treatment” group
denoted 1, and the ”Control” group, denoted 0.

#### Outcome

Outcome: Outcome is the recorded measure of HbA1c levels in mmol/mol 18-24
months after initiation of second line therapy.

#### Objectives

The estimand of interest is the Average Treatment Effect (ATE) of sustained
treatment with DPP4, compared to SU, over 18 months. Of secondary interest is
the how the the effect of treatment changes over 6 monthly intervals. Due to
the risk of unmeasured confounding, we use a measure of physicians
prescription preference (PP) taken over time as an IV. Full details are in the
methods section.

#### Covariates

Information on covariates measured at baseline is summarised in Table 3 in the
Appendix. Data is available on age, gender, ethnicity, HbA1c levels, Body Mass
Index, systolic blood pressure, blood lipid profiles, kidney function, and
history of stroke and hypoglycemic events. Co-prescription history of statins
and beta blockers was also available. No patients were prescribed Insulin in
this study. We consider that prescription preference may depend on some of
these confounders. Patients on DPP4 have lower HbA1c levels at baseline and
prior to 2nd line therapy, with higher levels of Body Mass Index over 34.
Notably, patients were majority non-white, with around 75% recorded as Black,
South Asian, or Other Ethnicity.

##  3 Methods

###  3.1 Overview

Suppose $T$ time periods for which we observe a treatment $A_{t}$ up to time
$T$. For the rest of the paper we will refer to ”treated” $(A_{t}=1)$ and
”control” $(A_{t}=0)$ groups. When referring to variables with no subscript
$t$, this is defined as the set of all observations of that variable, that is
$A=(A_{1},\ldots,A_{t})$ etc. We observe a continuous end of study outcome $Y$
at time $T+1$, and observed, and unobserved time varying confounders $L_{t}$
and $U_{t}$, confounding the effect of $A_{t}$ on $Y$. We define $Z_{t}$ as
the time-varying instrument.

We represent our expected data setups in the Directed Acyclic Graph (DAG)
shown in Figure 1. Previous works with time dependent variables [7, 6] assumed
similar data structures with the red directional arrows removed. We refer to
this as the ”Simple” data setup. A more ”Complex” setup includes these arrows,
and is akin to structures discussed in [5, 8].

Uniquely difficult aspects of this setup is the introduction of time-dependent
$Z-A$ confounding, direct associations between instrument times, and a much
more complicated $Z-A$ relationship, whereby $A_{t}$ is both directly
influenced by $Z_{t-1}$, and vice versa. These relationships may create
problems for standard time varying causal methods as we explain below.

Taking our motivating example, $Z_{t}$ is preference for prescribing DPP4 over
SU over the period $t$ (PP). We expect PP to be correlated over time. We might
also expect that certain time varying health related measures $L_{t}$,
progression of disease, blood glucose levels or patient self reported
satisfaction, could drive a change in which drug is preferred over time beyond
new emerging research and word of mouth. Moreover, PP is likely to be dictated
by prescriptions given over the long term, particularly for patients with a
long history of treatment. Arrows running from past PP to future treatment may
occur when repeat prescriptions are handed out, and as such preference at
initiation remains a direct factor on assignment. We therefore expect that two
stage methods to handle the complex setup of the DAG are needed.

$L_{1}$$L_{2}$$L_{3}$$Z_{1}$$A_{1}$$Z_{2}$$A_{2}$$Z_{3}$$A_{3}$$Y$$U_{1}$$U_{2}$$U_{3}$
Figure 1: DAG of Complex data setup with $T=3$ treatment periods

Define $Y(a)$ as the counterfactual outcome that would have been observed
under some treatment regime $a=(a_{1},\ldots,a_{T})$. A causal effect is
sought by estimating the contrast in counterfactual outcomes under different
treatment regimes. We define this contrast as the Average Treatment Effect
(ATE), the average difference in counterfactual outcome in the population when
always treated, versus never treated

$E[Y(1,1,1)-Y(0,0,0)].$

As with all causal methods, we make the assumptions of counterfactual
consistency and positivity throughout the paper [4].

In time fixed settings, a common alternative contrast is the Local Average
Treatment Effect (LATE). This is defined as the effect of treatment in those
for whom a change in IV status induces a change in treatment status. When $Z$
is defined as assignment to treatment, and $A$ treatment taken, this is
interpreted as the ATE in compliers, those who always take the treatment they
are assigned to. We do not consider the LATE in this paper, as it requires
certain monotonicity conditions on the data. These conditions, as well as what
exactly constitutes a complier, are not well defined, if at all, with a time-
varying treatment [23].

###  3.2 IV Assumptions

In time-varying settings we make the following specific assumptions for the
$Z_{t}$ to be valid instruments, and are multivariate extensions of those in
time fixed settings. Define $M_{t}$ as the set of variables for which, when
conditioned upon, $Z_{t}$ satisfies the following assumptions.

  1. 1.

IV1 (IV Relevance): There exists a measurement of $Z$ and $A$ at each time
point, so that we have as many instruments as treatment times, and there
exists an association between $Z_{t}$ and $A_{t}$, conditional on $M_{t}$, at
each time point. This association must sufficiently strong.

  2. 2.

IV2 (Exclusion Restriction): $Z_{t}\perp\\!\\!\\!\perp Y|M_{t}\;\forall t$.

  3. 3.

IV3 (Conditional Exchangeability): Provided there are no unmeasured
confounders between $Y$ and $Z$, we then have that $Z_{t}\perp\\!\\!\\!\perp
U|M_{t}\;\forall t$, that is $Z$ independent of any other unmeasured
confounding bewteen $A$ and $Y$. IV2 and IV3 can be formalised as the
assumption $Y(a)\perp\\!\\!\\!\perp Z_{t}|M_{t}\;\forall\;t\;and\;a$.

For our DAG we have in the simple case $M_{t}=Z_{t-1}$, and for the complex
case $M_{t}=(A_{t-1},Z_{t-1},L_{t})$. We define $M$, as the union of these
sets for all $t$.

###  3.3 Subtantive Models

Under these assumptions, established works in 2SLS settings present the
substantive model of interest as a Linear Structural Mean Model [17], which
models the observed outcome Y as

$E(Y|A,Z,M)=E[Y_{0}|A,Z,M]+\sum_{t=1}^{T}\beta_{t}A_{t}$ (1)

where $Y_{0}=Y(0,\ldots,0)$, and our interest lies in the effects of $A_{t}$
on outcome $\beta_{t}$. We cannot estimate the ATE under this model without
additional assumptions. Firstly, we make the assumption in this paper that the
treatments $A_{t}$ affect $Y$ linearly, with no interaction terms between
treatment times. Secondly, we require that the shift in conditional mean
outcome is the same in the treated arm, as in the control arm

$E[Y(\bar{a}_{t},0)-Y(\bar{a}_{t-1},0)|\bar{A}_{t-1}=\bar{a}_{t-1},A_{t}=1,\bar{Z}_{t},M_{t}]=E[Y(\bar{a}_{t},0)-Y(\bar{a}_{t-1},0)|\bar{A}_{t-1}=\bar{a}_{t-1},A_{t}=0,\bar{Z}_{t},M_{t}].$

for all $t$, where $\bar{A}_{t}$ represents the history of treatment up to
time t. This condition was described in [24] and more recently in [8]. It can
be shown that this property holds under the ”no current treatment interaction
property”, that is there is no effect modification of the effect of $A_{t}$ on
$Y$ by other variables, a property that clearly holds when the casual effect
is linear in $A$. When this property holds, the ATE can be estimated as the
sum of the coefficients $\beta_{t}$.

###  3.4 Two Stage Least Squares Methods

### Standard 2SLS

The standard 2SLS methodology that is well known in time fixed situations [1]
was generalised to the case of multiple exposures periods (or multiple
separate exposures) via Multivariate Mendelian Randomisation (MVMR) [15, 4].
This applied 2SLS to the case of a time varying exposure, with at least as
many measured IVs at baseline. In this paper we have a just-identified
situation, where we have as many instruments as time periods. Standard 2SLS
can be performed by the following steps.

  1. 1.

First Stage Models: For each $t$, postulate and fit a series of first stage
models for each $A_{t}$ $f_{A_{t}}=P(A_{t}|Z,L^{M})$ by fitting a main effects
Ordinary Least Squares (OLS) model for $A_{t}$ on all instruments $Z_{t}$ and
confounders for which $Z_{t}$ depend $L^{M}$.

  2. 2.

Obtain predicted values for each $A_{t}$ from the fitted models, defined as
$\hat{A_{t}}=E(A_{t}|Z_{1},\ldots,Z_{t},L^{M},\alpha^{*}_{t})$, where
$\alpha^{*}_{t}$ are the fitted coefficients.

  3. 3.

Second Stage Model: Postulate and fit a main effects OLS model for $Y$ against
the predicted values of $\hat{A}_{t}$ and $L^{M}$.

$f_{Y}=E(Y|\hat{A}_{1},\ldots,\hat{A}_{t};\beta)=\beta_{0}+\sum_{i=1}^{T}\beta_{t}\hat{A}_{t}+\beta_{L}L^{M}.$

using OLS regression. The estimated coefficients $\beta_{t}^{*}$ of the fitted
model are the causal effects of interest of Equation (1), for which the ATE
can be derived.

With time-varying IVs we expect that standard 2SLS can attain unbiased
estimates of causal effects in the simple data setup of Figure 1, without the
red arrows. In this instance, the $Z_{t}$, can be thought of functionally as
correlated IVs at baseline, which was demonstrated to perform well in recent
work [4]. We anticipate however that standard 2SLS may be unsuited to the more
complex data setup.

Firstly, we suspect that 2SLS may not be able to address time-varying
confounding between the $Z_{t}$, $A_{t}$ and $Y$. In cases of confounding only
at baseline, 2SLS can control for this by including any baseline confounders
in the first and second stage models [17]. With time-varying confounding,
conditioning on confounders beyond baseline blocks indirect causal effects of
$A$ on $Y$ through $L_{t}$, as well as risking collider bias. As such, the
third IV assumption, exchangeability cannot be satisfied as a uncontrolled for
pathway exists from $Z_{t}$ to $Y$ via $L_{t}$.

A second problem lies with the associations between $A$ and $Z$, as they
induce collider bias into the first stage models. In Figure 1 for example,
once the first stage model for $A_{1}$ conditions on $Z_{2}$, it collides, and
subsequently bias, the $Z_{1}-A_{1}$ association.

### Robust 2SLS (R2SLS) For Time Varying Instruments

To remedy this problem we now present the main method of the paper. We turn to
Robust 2SLS methodologies [25, 17, 18]. These methods deal with instances in
time fixed scenarios where exhangeability of $Z$ is dependent on some
unmeasured confounder at baseline $V$. These methods propose that before the
first and second stage model is fitted, a parametric model for the instrument
$Z$ is postulated

$f_{Z}=E(Z|V).$

A estimating equation method is then described that equates to the following
procedure.

  * •

Fit model $f_{Z}$, and calculate predictions of $Z$ as $Z^{*}$

  * •

Attain the residuals of $f_{Z}$, $Z-Z^{*}$

  * •

Perform standard 2SLS, with $Z-Z^{*}$ in place of $Z$, controlling for $V$ in
the first and second stage models.

To handle the complex data structure of Figure 1, we consider the following
extension of these methods to a time varying setting.

  1. 1.

Postulate a model for each $Z_{t}$ given $M_{t}$

$f_{Z_{t}}=E(Z_{t}|M_{t};\gamma)$

using an appropriate regression model. Fit this model and calculate
predictions for $Z_{t}$ as

$\hat{Z}_{t}=E(Z_{t}|M_{t};\gamma^{*})$

.

  2. 2.

From this define residuals $Z_{t}^{res}=Z_{t}-\hat{Z}_{t}$

  3. 3.

Now perform 2SLS, replacing $Z_{t}$ with $Z_{t}^{res}$.

  4. 4.

Standard 2SLS may be perform without conditioning on $L_{t}$ in the first and
second stage models.

One way to think of this method is that it is an application of standard 2SLS,
taking $Z_{t}^{res}$ as the time-varying IVs of the method. These residuals
instruments are sometimes known as the Frisch-Waugh-Lovell residualised
instruments, that are designed to allow the fitting of partial regression
models, and are orthogonal (unassociated) with $M_{t}$. This property is
crucial for two reasons. Firstly it ”models out” the associations between
$A_{t-1}$ and $Z_{t}$, and thus removes the collider bias they induce in the
first stage model.

Secondly we have that $Z^{res}_{t}|L_{t}$ and thus there is a choice as to
whether to control for $L^{M}$ in the first and second stage models at all. In
the complex case this means we no longer need to control for time-varying
confounding in the first and second stage models which blocked indirect
effects of $A_{t}$ on $Z_{t}$. We can thus obtain unbiased estimates of the
causal effects in the complex case provided that $f_{Z_{t}}$ is correctly
specified.

These methods have been referred to as g-estimators [18], and in fact R2SLS is
equivalent to to the g-estimation procedure recently considered in [4], and
[8] for time varying instruments, when the first and second stage are fit
using OLS models. We show this in the appendix.

### Robustness properties of R2SLS

We take specific note of the use of R2SLS in the simpler data setup. Recent
work [17] stress tested this methodology in the single time period case, in
situations in which $A_{t}$ and $Y$ depend on non linear relationships, such
as squared terms of $L_{t}$ and noted that the method was doubly robust.
Provided that either $f_{Z_{t}}$ or $f_{Y}$ was correctly specified.

This robustness however does not extend to the complex case, where we use
R2SLS as we cannot control for confounding beyond baseline in $f_{Y}$. In
essence we are leveraging the doubly robust property to avoid modeling
problematic associations in the first and second stage models. But in doing
so, we are reliant on correct specification of $f_{Z_{t}}$. How vulnerable
this method is to misspecification of $f_{Z_{t}}$ is a point of interest.

One consideration is what type of linear models should be used. We consider
situations in with $f_{Z_{t}}$ is fit by its natural linear model and the
first and second stage are fit by Ordinary Least Squares (OLS). When $A_{t}$
is binary, this corresponds to a Linear Probability Model (LPM). The appendix
considers other popular choices of first stage model in the case of binary
treatment.

##  4 Simulations

###  4.1 Data generating mechanism

We will test standard 2SLS and R2SLS via simulation. The major objectives of
the simulation are firstly, to determine the consistency and efficiency of
standard 2SLS, R2SLS-L and R2SLS under the data setups of Figure 1 under
varying IV strengths and sample sizes. Secondly, we aim to test the robustness
of R2SLS in simple and complex situations in which the $Z_{t}$, $A_{t}$ and
$Y$ depend on square terms of confounders, or $Z-A$ interaction terms, leading
to misspecification of $f_{Z_{t}}$. The data is simulated as follows.

  * •

The time-varying instrument $Z_{t}$ is simulated as
$logit(P(Z_{t}=1))=\mu_{Z_{t}}-E(\mu_{Z_{t}})$ where in the simple case

$\mu_{Z_{t}}=Z_{t-1}$

and the complex case as

$\mu_{Z_{t}}=Z_{t-1}+A_{t-1}+3*L_{t}+\sigma_{Z}(L_{t}^{2}).$

Here $\sigma_{Z}$ introduces typical terms that may trigger misspecification
of $f_{Z_{t}}$.

  * •

$A_{t}$ is a binary variable, generated using

$logit(P(A_{t}=1))=\Phi(\mu_{A_{t}}-E(\mu_{A_{t}}))(1-\Delta_{t})+Z_{t}\Delta_{t}.$

where

$\mu_{A_{t}}=A_{t-1}+L_{t}+U_{t}$

in the simple case and

$\mu_{A_{t}}=Z_{t-1}+A_{t-1}+L_{t}+U_{t}+\sigma_{A}(L_{t}^{2})$

in the complex case with $\Phi$ denoting the standard normal cdf. We use
$\sigma_{A}$ to introduce non linear terms and interactions into the models.
We simulate this way based on work in [6, 8] where $\Delta_{t}$ can be chosen
to control the strength of the instrument $Z_{t}$, simulated as
$\Delta_{t}=\Phi(\alpha)$.

  * •

Lastly $Y$ is normally distributed with variance 1, and mean

$E(Y)=\sum_{i=1}^{T}(U_{t}+A_{t}+L_{t})+\sigma_{Y}L_{1}^{2}$

The true values for $\beta_{t}$ are $(3,2,1)$, and subsequently the ATEs true
value is 6.

###  4.2 Implementation

In each simulation we generate 1000 datasets. Confidence intervals (CIs) are
obtained via a percentile bootstrap method using 1000 bootstrapped datasets.
In all simulations standard 2SLS is performed as in section 3, controlling for
$L_{1}$. R2SLS is performed with $f_{Z_{t}}$ defined as main effects logistic
regression model on $Z_{t-1}$ for the simple case and, in the more complex
case $Z_{t-1}$, $A_{t-1}$ and $L_{t}$. $f_{A_{t}}$. The first stage model is a
main effects LPM regression of $A_{t}$ on $Z^{res}$, and the second stage
model is a main effects OLS regression of $Y$ on $\hat{A}$

We look to vary the sample size $n$ which will be set to 1000 or 5000, and the
parameter $\alpha$, which can be set to values between 0 and 1 to influence
the strength of association between $Z_{t}$ and $A_{t}$. We choose these
values as $0.1$, $0.3$ and $0.5$, which correspond to weak, moderate, and
strong associations respectively.

We also look to test the robustness of R2SLS in the complex case, which cannot
use doubly robust methods to misspecified models, when we set $\sigma_{Z}$ and
$\sigma_{Y}$ to $0$ or $1$ respectively. When non-zero, these terms introduce
non linear terms into $f_{Z}$ and $f_{Y}$, specifically square terms of the
baseline confounder $L_{1}$ and interactions between past treatment and
instrument, leading to misspecified models.

As recommended in [26] we report for the ATE the average absolute bias, Root
Mean Square Error (RMSE) and Coverage defined as the proportion of datasets in
which the bootstrap CI included the true value.

###  4.3 Results

Tables 1 and 2 present the results of the simulations in the case of no non-
linear terms $(\sigma_{Z},\sigma_{A},\sigma_{Y})=(0,0,0)$. The full simulation
results are in the Appendix, Tables 4, 5 and 6. We can see that standard 2SLS
shows the same effectiveness as the robust case in the setup for our simpler
simulation. This is expected, and identified by the work in [4]. With the
simple setup of Figure 1, the $Z_{t}$ can be treated as though they were
baseline instruments. Results show low bias and good coverage in all cases
where instrument strength is acceptably strong. In table 4 2SLS was employed
using a first stage probit model rather that OLS, which also performed well.
Results began to show bias only at low correlation strengths $\alpha=0.1$ and
low sample size $n=1000$.

As expected in complex scenarios with time varying confounding R2SLS is the
only method capable of obtaining unbiased estimates conditional on sufficient
instrument strength. 2SLS shows high bias in all cases. Probit based first
stage models in Table 5 no longer performed well, even for R2SLS. This may
infer that the orthogonality to $M_{t}$, for which R2SLS is based upon is not
guaranteed when using a non OLS model. It also suggests that unaccounted for
time varying confounding can lead to misleading results. We can infer from
Tables 1 that R2SLS with LPMs is the only reliable method available in complex
circumstances.

Table 2 shows the results using R2SLS for the complex data setup when we add
non linear terms to the models. We add terms into the models for $Z$ and $A$
at the same time, in order to control the instrument strength. As expected,
when all 3 models are misspecified we encounter biases and less efficient
results in all cases. This bias is quite low however in cases of strong
instrument strength and sample size.

When $\sigma_{Z}=1$ and $\sigma_{Y}=0$ results show bias similar to when both
are 1. This is as expected as the inability to model the Y-M association means
there is not robust protection from misspecified $f_{Z_{t}}$. The opposite way
around however results are unbiased, subject to instrument strength, as
correct specification of $f_{Z_{t}}$ protects from misspecification of
$f_{Y}$, as expected.

n | $\alpha$ | Bias | RMSE | MCE | Coverage
---|---|---|---|---|---
Simple Simulation
Standard 2SLS
5000 | 0.5 | 0.004 | 0.42 | 0.006 | 93.5
5000 | 0.3 | 0.012 | 0.52 | 0.009 | 93.6
5000 | 0.1 | 0.096 | 0.85 | 0.022 | 94.2
1000 | 0.5 | 0.013 | 0.63 | 0.012 | 94.5
1000 | 0.3 | 0.048 | 0.79 | 0.020 | 95.6
1000 | 0.1 | 0.128 | 4.18 | 0.554 | 98.9
R2SLS (g-estimation)
5000 | 0.5 | 0.004 | 0.42 | 0.006 | 93.5
5000 | 0.3 | 0.012 | 0.52 | 0.009 | 93.6
5000 | 0.1 | 0.096 | 0.85 | 0.022 | 94.2
1000 | 0.5 | 0.013 | 0.63 | 0.012 | 94.5
1000 | 0.3 | 0.048 | 0.79 | 0.020 | 95.6
1000 | 0.1 | 0.128 | 4.18 | 0.554 | 98.9
Complex Simulation
Standard 2SLS
5000 | 0.5 | 1.900 | 1.38 | 0.003 | 0.0
5000 | 0.3 | 2.513 | 2.81 | 0.238 | 76.3
5000 | 0.1 | 3.397 | 1.91 | 0.041 | 2.7
1000 | 0.5 | 1.895 | 1.38 | 0.007 | 0.2
1000 | 0.3 | 2.781 | 4.12 | 0.531 | 78.3
1000 | 0.1 | 3.361 | 2.87 | 0.240 | 40.9
R2SLS (g-estimation)
5000 | 0.5 | 0.010 | 0.33 | 0.003 | 96.1
5000 | 0.3 | 0.014 | 0.42 | 0.006 | 95.2
5000 | 0.1 | 0.638 | 3.21 | 0.325 | 96.2
1000 | 0.5 | 0.006 | 0.50 | 0.008 | 95.8
1000 | 0.3 | 0.023 | 0.644 | 0.013 | 95.4
1000 | 0.1 | 5.352 | 11.439 | 4.140 | 98.7
Table 1: Simulation results, targeting the ATE with coverage based on b=1000 bootstrapped samples. $\sigma_{Z}$, $\sigma_{A}$ | $\sigma_{Y}$ | n | $\alpha$ | Bias | RMSE | MCE | Coverage
---|---|---|---|---|---|---|---
R2SLS (g-estimation
1 | 1 | 5000 | 0.5 | 0.026 | 0.60 | 0.011 | 94.2
1 | 1 | 5000 | 0.3 | 0.083 | 0.75 | 0.017 | 93.7
1 | 1 | 5000 | 0.1 | 0.690 | 2.22 | 0.154 | 98.3
1 | 1 | 1000 | 0.5 | 0.101 | 0.90 | 0.026 | 95.5
1 | 1 | 1000 | 0.3 | 0.192 | 1.17 | 0.042 | 94.6
1 | 1 | 1000 | 0.1 | 6.922 | 12.41 | 4.868 | 98.3
1 | 0 | 5000 | 0.5 | 0.013 | 0.63 | 0.012 | 93.7
1 | 0 | 5000 | 0.3 | 0.016 | 0.77 | 0.018 | 94.6
1 | 0 | 5000 | 0.1 | 0.593 | 2.36 | 0.176 | 97.6
1 | 0 | 1000 | 0.5 | 0.075 | 0.95 | 0.028 | 94.0
1 | 0 | 1000 | 0.3 | 0.165 | 1.22 | 0.047 | 94.6
1 | 0 | 1000 | 0.1 | 7.000 | 12.42 | 4.877 | 98.2
0 | 1 | 5000 | 0.5 | 0.010 | 0.33 | 0.003 | 95.7
0 | 1 | 5000 | 0.3 | 0.015 | 0.42 | 0.006 | 95.1
0 | 1 | 5000 | 0.1 | 0.228 | 1.45 | 0.066 | 96.5
0 | 1 | 1000 | 0.5 | 0.002 | 0.50 | 0.008 | 95.7
0 | 1 | 1000 | 0.3 | 0.01 | 0.64 | 0.012 | 95.5
0 | 1 | 1000 | 0.1 | 3.440 | 12.20 | 4.709 | 98.8
Table 2: Simulation results targeting the ATE, using R2SLS (g-estimation) for
complex data setups with mispecified models for $Z$ $A$ and $Y$.

##  5 Case Study

###  5.1 Instrument Definition: GP Prescription Preference

Our instrument is a measure of GP prescription preference (PP) for DPP4 over
SU over time. There are 139 GPS in the data, with an average of 27 patients
each, ranging from 1 to 98 patients. A recent paper [11] summarized well the
various measures to approximate GP preference via proportion of prescriptions
issued. Available prescription data does not include specific dates, so
subject specific measures of PP are not calculable. We instead consider GP
specific measures of preference at each 6 month period based on the
definitions of $IV_{prop}$ [11].

  1. 1.

$PP^{Cal}$: A measure of GP preference at each 6 monthly calendar period is
taken as the proportion of all prescriptions of DPP4 as second line treatment
within that 6 monthly calendar period. An individuals value of $PP^{Cal}$ at
some follow up time $t$ is then GP preference during the 6 months calendar
period when $t$ occured. We label this $PP_{t}^{Cal}$

  2. 2.

$PP^{t}$: Alternatively, a measure of PP over follow up time $t$, rather than
calendar time can be considered. The proportion of a GPs prescriptions at
initiation are taken as $PP_{1}^{t}$. This is repeated for follow up times 2
and 3. This represents a measure of how GP may change preference based on how
long a patient has been taking second line treatment.

A visual representation of the definitions of $PP$ is shown found in Figure 3
in the Appendix. Preference taken at 6 month period shows a clear shift in
preference over time towards DPP4, with a proportion of around 15% in 2013 to
40% in 2017 . However, $PP^{t}$ shows minimal change between follow up
periods.

#### IV1: IV Relevance

It is important in a practical context to investigate if $Z_{t}$ satisfy the
three main IV assumptions. Unfortunately, only IV1, IV relevance can be tested
from the available data. Few patients switch from one treatment to another
over the study period. We need to investigate if there is firstly sufficient
strength of association between PP and treatment at initiation, and then
assess if there is strong enough change in this association over time for
multiple instruments.

In a time fixed setting, the Cragg-Donaldson F-statistic is used to evaluate
the strength of instruments [27, 28, 29]. An F-statistic of 10 is historically
cited as a sign of a sufficiently strong instrument, however more recent work
suggesting a value of 100 may be necessary [30]. This test is not however
appropriate to evaluate the overall strength of multiple IVs. The Sanderson-
Windmaijer conditional F-test [15, 31], adapts the Cragg-Donaldson statistic
to provide a value for each $A_{j}$, testing the remaining predictive power in
the $Z_{j}$ once all other $A_{j}^{res}$ have been calculated in the first
stage model. Analysis of standard errors, unusual estimates, and inspecting
the influence of $Z_{t}$ in the first stage models, can also infer whether
there is sufficient instrument strength.

The first stage F-statistics for PP at $t=1$ were 1422 and 927 for $PP^{Cal}$
and $PP^{t}$ respectively. Strengths similar to this were found in a recent
paper on PP in assessment of T2D second line treatment [32]. The conditional
F-tests for $PP^{Cal}$ and $PP^{t}$ however ranged between 30 and 70, below
the recently recommended standard of 100, and dramatically lower then the
first stage F-statistics. Correlations between PP times were very high, above
0.85 for $PP^{Cal}$, and in the case of $PP^{t}$ they neared 1. This explains
the heavy drop in predictive power and indicates a high risk of collinearity.

#### IV2: Exclusion Restriction

It cannot be determined from the available data if there is a direct
association between GP preference and HbA1c levels at any future time. One
possible avenue is that GPs with a preference for SU over DPP4, may be of
”higher standard”, leading to an improvement in HbA1c levels. Whilst this
assumption cannot be tested, it can be reasoned that standard of care between
practices will depend on numerous other more critical factors such a funding,
region and standards of follow up care, and such a possibility is unlikely.

#### IV:3 Exchangeability

GP preference may be affected by past confounders, and measurements of
outcome. Dependence on past preference and treatment history is likely, but
can be easily controlled. However a GP identifying poor past performance of
patients on one drug may switch preference as a result. This opens the
possibility of preference drifting over time towards DPP4 being dependent on
other health related measures other than past treatment, such as HbA1c
history. We do not for example, have full information on comorbidities.
However, the preference of a GP may be reinforced by the overall health of
their previous patients prior to that 6 month period. We argue that health
related issues are normally considered on a patient by patient basis by GPs,
and thus should directly affect treatment assignment, rather directly
impacting preference. It is also unlikely for exchangeability to be met at
initiation, only to fail at later follow up times, as health related factors
are always being considered. It is most crucial therefore to be confident of
exhangeability at the first time period.

###  5.2 Estimating Approaches

We use 2SLS and R2SLS to investigate the effect of sustained treatment with
DPP4 versus SU on HbA1c levels at 18 months, estimating the Average Treatment
Effect and as a secondary measure, the trajectory of treatment effect over
time. We repeat these analyses for $PP^{Cal}$ and $PP^{t}$. The models of
interest are as follows.

  * •

Standard 2SLS: The Standard 2SLS as described in Section 3.4. This represents
a benchmark analysis, as we anticipate that it will attain biased causal
effects.

  * •

R2SLS: The Robust 2SLS methodology described in Section 3.4.

  * •

Ridge R2SLS: The R2SLS method in which the second stage model is fit via Ridge
regression methods to address collinearity. See Appendix for details

Confidence intervals are obtained by non parametric bootstrapped with $b=1000$
resampled datasets.

R2SLS models $Z_{t}$ against history of treatment, and instrument as well as
history of HbA1c levels, Median HbA1c levels prior to initiation, smoking
status, calendar period and ethnicity. These were identified as potentially
unbalanced by PP at $t=1$ (see Appendix, Table 8).

Ridge R2SLS applies ridge regularisation penalty for the second stage in an
attempt to deal with multicollinearity in the second stage model. An algorithm
is described in the appendix, based on ridge trace plot stabilising efforts,
to identify a penalty for which the $\beta_{j}$ were all within 0.01 of the
values at the previous tested penalty.

###  5.3 Results

A full table of results are shown in the Appendix, Table 7. We show the
results for the ATE in Figure 2. All methods suggest a reduction in HbA1c
levels with sustained treatment with DPP4 compared to SU. Standard 2SLS
identified a decrease in HbA1c levels at 18 months of around 5 mmol/mol. Using
R2SLS, this decrease was between 2.5 to 3 points. All methods indicated a
significant reduction in HbA1c levels.

There was a modest gain in precision with RS2LS over standard 2SLS, and a
reduction in the effect size. This may indicate that 2SLS was affected by
collider bias or unmeasured confounders of the IV. Results with $PP^{Cal}$
were also slightly more precise than with $PP^{t}$.

Whilst R2SLS was able to identify an stable estimate of the ATE, trajectories
however suffered from very unstable estimates and confidence intervals.
Inspection of the first stage models indicate that PP at time 1, explains the
vast majority of the treatment-PP association over time, leading to jointly
weak instruments. The second stage model also displayed very high
collinearity, with a Variance Inflation Factor (VIF) in the 100s or 1000s,
which will heavily inflate standard errors.

Use of Ridge R2SLS soothes the high instability caused by collinearity and
pulls estimates towards the null. The ATEs are similar with a modest reduction
in the width of the confidence interval. Trajectories imply a dampening of the
effect of DPP4 over time. However trajectories using this method are likely to
remain biased due to weakly joint instruments, an issue which Ridge
regularizations are unsuited to address. Stabilising estimates with Ridge
regularisation requires further work before being recommended more generally.
This is considered in the discussion.

Figure 2: Forest Plot of case study results for PP

##  6 Discussion

In this paper, we investigate the use of Standard and Robust 2SLS
methodologies in a full time-varying setting. We find by simulation that
Standard 2SLS can perform well in situations in which the IV has simple
relationships between treatments and confounders, but showed biased estimates
in more complicated scenarios. Robust 2SLS was found to perform well even in
complicated setups. Results for R2SLS showed bias only at weak correlation
strengths and low sample size. R2SLS also performed relatively well in when
non linear terms were introduced. An evaluation of the effect of sustained
treatment intensification with DPP4 versus SU over 18 months found a
significant reduction in HbA1c levels of between 2.5-5 mmol/mol.

This paper contributes to the literature in the following ways. We critically
evaluate the performance of 2SLS methodology in a full time-varying setting
with time-varying confounders, instruments, and complex relationships over
time. We adapt literature on robust 2SLS methods in time fixed situations to a
multivariate setting and demonstrate that this adaptation can successfully
recover time varying causal effects by modeling dependence on prior time
periods. We identify that R2SLS can be related to recent advances in
g-estimation for time-dependent instruments. We also assess the performance of
R2SLS in multivariate setting under non linear relationships between the IV,
treatment, and time varying confounders.

We build on prior literature focusing on Multivariate Mendelian Randomisation
[15] by extending its methodologies to time-varying instruments and time
varying-confounding. We consider practical situations rarely encountered in
genetic studies and find the methods can handle these data setups. We extend
the existing work on Robust 2SLS methods detailed for time fixed scenarios
[18, 17], and find that they are suitable to handle time varying confounding,
and retain the robust properties, identified in the univariate setting. We
also confirmed that the relative robustness to non linear relationships appear
to extend to multiple time periods, provided the model for the IV is correctly
specified.

Recent studies have compared SU versus DPP4 in both a clinical trials setting
[33, 34], and a retrospective cohort [35]. In particular, recent work using
IVs compared glycemic control in SU and DPP4 using physician preference as a
baseline instrument [32, 11]. In these studies, around 75% of patents were
white. In this paper, our study was 75% non-white, specifically patients of
South Asian and black ethnicity. Prevalence of diabetes is estimated to be two
to four times higher in Asian and Black Ethnic Groups [36], and there is as
such considerable interest in the external validity of previous studies on
glycemic control to non-white populations in the UK. We identify effects of
DPP4 versus SU that are in line with these previous studies, providing
external validity to these papers and evidence of the effectiveness of DPP4 in
non white populations.

Simulations showed that whilst non linear and interaction terms in model for
the instrument are relatively robust to bias there is, as is typical for IV
methods, a heavy penalty if the main effect of a confounder of the $Z-Y$
relationship is not controlled for. Unlike time fixed R2SLS the outcome model
can rarely be used as a means to make results doubly robust. Heteroskedastic
errors in the instrument may also impact the ability of the method to ”model
out” complex associations between $Z$ and $A$.

For the case study, we found that even with identifying PP with strong
strength at baseline, it was challenging to attain time-varying PP with enough
variation in association over time to identify treatment trajectories. This
was due to high correlations between the instruments. Specifically, PP at
initiation contained the majority of the predictive power of treatment
assignment in the future. Instruments showed a weakly sufficient conditional
F-statistic, however autocorrelated errors may have skewed this. Certain
confounding information was not available, and so we could only partly
investigate balance across the instrument. We also could not empirically test
the exclusion criteria.

Unstable estimates of the trajectory and highly collinear variables in the
second stage model were the result, which were partly addressed using Ridge
regression. However regularisation is unlikely to address the underlying
problem of weak instruments, and more sophisticated techniques are needed
going forward.

Further research may consider the use of techniques such as generalised least
squares to account for heteroskedastic error terms in the models for the IV
[37]. Difficulties with highly correlated IVs could be attacked in one of two
ways. Firstly, significantly improve the strength of the IV at time 1.
Secondly, minimise the correlations between the instruments. Near far matching
methodologies could identify pairs of GPs with the furthest possible
difference in preference at initiation, or who are most different between two
instrument times. We may also consider methods such as the Abrahamowics
method, [32] to attempt to model the point in which a physicians preference
changes, to construct multiple IVs.

In a prior work [8] prescription preference for Rheumatoid Arthritis treatment
was taken over time to some success. However patents switched treatment
strategies frequently, with a choice between more than two medications and the
conditional F-statistic was dramatically higher than any identified for this
study. This could indicate that time-varying IVs are more suited to studies
with more than two treatment options .

To conclude, we have demonstrated and provided simple guidance to use R2SLS in
a time varying instrument setting. The proposed method works well and is
readily accessible to users in a variety of fields. Implementing a time-
varying IV in practice is challenging due to high correlations between the IV
over time. This makes it difficult to find sufficiently strong instrument
strength which varies over time to estimate trajectories of causal effects.
The viability of PP as a time varying instrument remains one of important
research.

### Acknowledgments

The authors thank the Queen Mary University of London, Clinical Effectiveness
Group and Professor Bianca De Stavola for the data.

### Conflicts of Interest

There are no conflicts of interest to declare.

## Appendix

### Baseline Characteristics

Characteristics |  Sulfonylureas (n=1733) |  DPP4 (n=828)
---|---|---
Age (SD) | 56.8 (12.0) | 56.9 (11.7)
Male (%) | 53.8 | 48.5
Ethnicity (%) |  |
White | 25.3 | 24.9
South Asian | 50.4 | 53.7
Black | 19.1 | 16.5
Other | 5.1 | 4.8
Non-Smoker (%) | 28.0 | 27.7
HbA1c level (SD) (mmol/mol) | 67.7 (17.6) | 65.8 (15.5)
Median HbA1c Prior (SD) to Initiation (mmol/mol) | 64.2 (17.4) | 61.5 (15.4)
Body Mass Index (%) |  |
$<25$ | 19.1 | 17.3
$25-29$ | 38.3 | 29.9
$30-34$ | 25.4 | 25.5
$>34$ | 17.2 | 27.1
Years of Diabetes |  |
$<3$ yrs | 29.4 | 28.3
$3-5$ yrs | 20.8 | 18.1
$>5$ yrs | 49.7 | 53.5
Systolic Blood Pressure |  |
$<100$ | 1.61 | 0.60
$100-139$ | 83.8 | 82.8
$>140$ | 14.5 | 16.5
Total Serum Cholesterol (mmol/L) (SD) | 4.12 (1.05) | 4.14 (1.10)
Glomerular Filtration Rate (eGFR) |  |
$>60$ | 91.6 | 92.1
$30-59$ | 7.2 | 6.8
$<30$ | 1.0 | 0.9
No of Hypoglycemic events |  |
$0$ | 98.1 | 98.1
$1+$ | 1.9 | 1.9
Number of strokes |  |
$0$ | 98.9 | 99.2
$1$ | 1.1 | 0.8
Medications (%) |  |
Statins | 86.3 | 89.1
Beta-Blockers | 16.6 | 17.0
Table 3: Patient baseline characteristics at initiation of second line
therapy.

### R2SLS Relation to g-estimation

In g-estimation the substantive model of interest is a Structural Nested Mean
Model (SNMM)

$E[Y(\bar{a}_{t},0)-Y(\bar{a}_{t-1},0)|\bar{A}_{t}=\bar{a}_{t},\bar{Z}_{t},M_{t}]=\beta_{t}a_{t}\qquad
t=1,\ldots,T-1.$ (2)

where $Y(\bar{a}_{t},0)$ is the counterfactual outcome under a patient’s
observed treatment history up to time $t$, and on control afterwards.

Under the assumptions of no current treatment interaction, and a linear
relationship bewteen $A$ and Y (see [8]) the SNMM can be reduced to a Marginal
Structural Model of the form

$E[Y(a)]=\alpha+\sum_{t=1}^{T}\beta_{t}a_{t}.$

Notably this model implies the Linear Structural Mean Model of 2SLS, which
models the observed outcome Y as

$E(Y|A,Z,M)=E[Y_{0}|A,Z,V]+\sum_{t=1}^{T}\beta_{t}A_{t}$

In order to prove the equivalence of RS2LS and G-estimation with time varying
IVS, we demonstrate the equivalence of their closed form solutions.

Define

$A=\begin{bmatrix}A_{11}&\ldots&A_{T1}\\\ \vdots&\vdots&\vdots\\\
A_{1n}&\ldots&A_{Tn}\\\ \end{bmatrix}$

and define $Y$, $\hat{A}$, $Z$ equivalently. For notational convenience we
label $Z^{res}=Z$ G estimation has the closed form solution.

$\beta^{gest}=Y^{{}^{\prime}}ZZ^{{}^{\prime}}A(A^{{}^{\prime}}ZZ^{{}^{\prime}}A)^{-1}$

and the robust 2SLS method can be shown to have the closed form solution

$\beta^{2SLS}=(A^{{}^{\prime}}PA)^{-1}(A^{{}^{\prime}}PY)$

where P is the projection matrix

$P=Z(Z^{{}^{\prime}}Z)^{-1}Z^{{}^{\prime}}$

and thus the transpose can be written as

$\beta^{2SLS^{\prime}}=[(A^{{}^{\prime}}PA)^{-1}(A^{{}^{\prime}}PY)]^{{}^{\prime}}$
$=(A^{{}^{\prime}}PY)^{{}^{\prime}}[(A^{{}^{\prime}}PA)^{-1}]^{{}^{\prime}}$
$=Y^{{}^{\prime}}P^{{}^{\prime}}A[(A^{{}^{\prime}}PA)^{{}^{\prime}}]^{-1}$
$=Y^{{}^{\prime}}PA[(A^{{}^{\prime}}PA)]^{-1}$

Due to $P$ being symmetric.

We see that the two closed form solutions are similar with $P$ replacing
$ZZ^{{}^{\prime}}$, and are in fact equal as we have, expanding $P$

$\beta^{2SLS^{\prime}}=Y^{{}^{\prime}}Z(Z^{{}^{\prime}}Z)^{-1}Z^{{}^{\prime}}A[(A^{{}^{\prime}}Z(Z^{{}^{\prime}}Z)^{-1}Z^{{}^{\prime}}A)]^{-1}$

noting that $Z^{{}^{\prime}}A$, $(Z^{{}^{\prime}}Z)^{-1}$ and
$A^{{}^{\prime}}Z$ are square matrices, we can rewrite the inverse term as

$=Y^{{}^{\prime}}Z(Z^{{}^{\prime}}Z)^{-1}Z^{{}^{\prime}}A(Z^{{}^{\prime}}A)^{-1}(Z^{{}^{\prime}}Z)(A^{{}^{\prime}}Z)^{-1}$
$=Y^{{}^{\prime}}Z(A^{{}^{\prime}}Z)^{-1}$
$=Y^{{}^{\prime}}Z(Z^{{}^{\prime}}A)(Z^{{}^{\prime}}A)^{-1}(A^{{}^{\prime}}Z)^{-1}$
$=Y^{{}^{\prime}}ZZ^{{}^{\prime}}A(A^{{}^{\prime}}ZZ^{{}^{\prime}}A)^{-1}=\beta^{gest}$

Thus these two methods are equivalent, under the assumption that there are as
many $Z$ as $A$, both the first and second stage models are fit by OLS, and
that $Z^{{}^{\prime}}A$, $(Z^{{}^{\prime}}Z)^{-1}$ and $A^{{}^{\prime}}Z$ are
invertable matrices.

### Penalised Regression Methods for Second Stage Model

One challenge that may be encountered with time dependent instruments is that
of collinearity. In situations in which a treatment does not vary much over
time, the $A_{t}$ will be highly correlated. If the $Z_{t}$ also show only
minor change over time, we are faced with the possibility that the majority of
the predictive power of the $Z_{t}$ lies within the first time period. That is
to say that each $A_{t}$ is far more strongly associated with $Z_{1}$ than any
other $Z_{t}$.

Under these circumstances, the predicted $\hat{A}_{t}$ may be highly
correlated and the second stage model may suffer from multicollinearity. This
can lead to highly unstable estimates of the time specific causal effects
$\beta_{t}$ and large standard errors.

To combat this we consider fitting the second stage model using Ridge
regression methods [38, 39]. Ridge regression fits the OLS model with penalty
function determined by a tuning parameter $\lambda$. The solution to the
second stage model is then.

$\beta=min_{\beta}\\{\sum_{i=1}^{n}(Y-\alpha-\beta_{t}\hat{A_{t}})^{2}+\sum_{t=1}^{T-1}\lambda\beta_{t}^{2}\\}$
(3)

which is better known in its closed form solution

$\beta=(\hat{A}^{t}\hat{A}+\lambda I)^{-1}\hat{A}^{t}Y$

with $I$ the identity matrix of appropriate dimension. We see that the penalty
$\lambda$ penalises large values of $\beta_{t}$ reducing instability in the
estimates, and breaking the collinearity of the $\hat{A}_{t}$. The compromise
is that it introduces bias to $\beta_{t}$. The larger the value of $\lambda$
the greater penalty to larger estimates of $\beta$, which biases estimates
towards zero. Thus the choice of $\lambda$ is important. To minimise bias, we
seek to choose the smallest $\lambda$ such that the colinearity in the second
stage model is suitably broken to provide stable estimates of $\beta_{t}$.
This is often done using ridge trace plots to observe the range in which
estimates stabilise but can be formalised with the following algorithm.

  1. 1.

Select a suitably wide range of values of $\lambda$, and take a set of $k$
increasing values $\lambda_{k}$ over this range

  2. 2.

Perform 2SLS using Ridge regression in the second stage with each penalty
$\lambda_{k}$. For each $k$, calculate the change in each parameter from the
previous run. $\|\beta^{k}-\beta^{k-1}\|$

  3. 3.

Select the first parameter $\beta^{k}$ for which this difference is within
some chosen tolerance $tol$ for each $\beta_{t}$

The tolerance should be small enough to attain stable estimates, but too small
will lead to choosing larger penalties and shrink estimates. Penalised
regression methods have been considered in the first stage model in the
context of Multivariate Mendelian randomisation to select from a high number
of genetic instruments [40], but to our knowledge not in this context in the
second stage model. We favour Ridge over other types of penalised regression
because it cannot shrink estimates to 0.

### Additional Simulations

#### Probit based First Stage Models for Binary Data

As with most theory on 2SLS, we have made the assumption that the first and
second stage models are fit using Ordinary Least Squares (OLS). This creates a
dilemma when the treatment and/or outcome are binary. Our approach is to fit
an OLS model to the binary data, known as a linear probability model (LPM).
LPMs however may not fit such data particularly well, where we might be expect
the distribution near probabilities close to 0 or 1 to become non linear.

One commonly explored avenue is to fit the first stage models using a probit
or logistic regression model. These are often referred to as Two Stage
Predictor Substitution (2SPS) [41] or ”forbidden” regression models [42].
Whilst non OLS models may fit the binary data better, their consistency is
still a matter of debate.

In the case of a binary treatment, non OLS first stage models cannot guarantee
that the residuals of the model are uncorrelated with the predicted values
$\hat{A}_{t}$. One noted trick however [42] is as follows

  1. 1.

Use 2SLS to obtain predictions of the treatment $\hat{A}_{t}$ from a first
stage model with chosen non OLS link function.

  2. 2.

Now perform 2SLS with an OLS first and second stage model, substituting the
instruments $Z_{t}$ with the predictions $\hat{A}_{t}$

This has been shown to be an efficient and consistent alternative to other
2SLS approaches [42].

In this paper we will only consider continuous outcomes. Two stage methods
with a binary outcome are a greater challenge. Methods using a second stage
LPM have been shown to be mainly unbiased, but most non LPM methods appear to
be inconsistent [42, 41], and depends on rarity of the outcome. This is beyond
the current scope of the paper. We repeat the simulations of the main body of
the paper, whilst additionally testing R2SLS in the case of a probit first
stage model, and using the ”trick” method described above. In the simpler
scenario, results show good or even slightly better performance that 2SLS or
R2SLS with an OLS model, but show notable bias in the complex case. This
indicates that the residualisation of the instruments is not effective when
not applying OLS models.

n | $\alpha$ | Bias | RMSE | MCE | Coverage
---|---|---|---|---|---
Standard 2SLS
5000 | 0.5 | 0.004 | 0.42 | 0.006 | 93.5
5000 | 0.3 | 0.012 | 0.52 | 0.009 | 93.6
5000 | 0.1 | 0.096 | 0.85 | 0.022 | 94.2
1000 | 0.5 | 0.013 | 0.63 | 0.012 | 94.5
1000 | 0.3 | 0.048 | 0.79 | 0.020 | 95.6
1000 | 0.1 | 0.128 | 4.18 | 0.554 | 98.9
R2SLS
5000 | 0.5 | 0.004 | 0.42 | 0.006 | 93.5
5000 | 0.3 | 0.012 | 0.52 | 0.009 | 93.6
5000 | 0.1 | 0.096 | 0.85 | 0.022 | 94.2
1000 | 0.5 | 0.013 | 0.63 | 0.012 | 94.5
1000 | 0.3 | 0.048 | 0.79 | 0.020 | 95.6
1000 | 0.1 | 0.128 | 4.18 | 0.554 | 98.9
R2SLS (Probit)
5000 | 0.5 | 0.027 | 0.42 | 0.006 | 94.0
5000 | 0.3 | 0.022 | 0.52 | 0.009 | 93.5
5000 | 0.1 | 0.054 | 0.84 | 0.022 | 94.1
1000 | 0.5 | 0.032 | 0.62 | 0.012 | 94.4
1000 | 0.3 | 0.040 | 0.78 | 0.020 | 94.4
1000 | 0.1 | 0.400 | 1.96 | 0.121 | 98.9
R2SLS (Probit ”Trick”)
5000 | 0.5 | 0.003 | 0.42 | 0.006 | 93.5
5000 | 0.3 | 0.008 | 0.52 | 0.009 | 93.6
5000 | 0.1 | 0.049 | 0.84 | 0.022 | 94.1
1000 | 0.5 | 0.006 | 0.63 | 0.012 | 94.2
1000 | 0.3 | 0.023 | 0.78 | 0.020 | 95.1
1000 | 0.1 | 0.342 | 1.56 | 0.076 | 99.9
Table 4: Simulation results for Simple data setup of Figure 1, based on b=1000 bootstrapped samples. LPM and Probit indicate first stage models fit using OLS or a probit regression model respectively. Probit ”Trick” denotes the method in which first stage predictions of $A_{t}$ using probit are subsequently used as instruments in 2SLS.  n | $\alpha$ | Bias | RMSE | MCE | Coverage
---|---|---|---|---|---
Standard 2SLS
5000 | 0.5 | 1.900 | 1.38 | 0.003 | 0.0
5000 | 0.3 | 2.513 | 2.81 | 0.238 | 76.3
5000 | 0.1 | 3.397 | 1.91 | 0.041 | 2.7
1000 | 0.5 | 1.895 | 1.38 | 0.007 | 0.2
1000 | 0.3 | 2.781 | 4.12 | 0.531 | 78.3
1000 | 0.1 | 3.361 | 2.87 | 0.240 | 40.9
R2SLS
5000 | 0.5 | 0.010 | 0.33 | 0.003 | 96.1
5000 | 0.3 | 0.014 | 0.42 | 0.006 | 95.2
5000 | 0.1 | 0.638 | 3.21 | 0.325 | 96.2
1000 | 0.5 | 0.006 | 0.50 | 0.008 | 95.8
1000 | 0.3 | 0.023 | 0.644 | 0.013 | 95.4
1000 | 0.1 | 5.352 | 11.439 | 4.140 | 98.7
R2SLS (Probit)
5000 | 0.5 | 0.535 | 0.74 | 0.004 | 1.2
5000 | 0.3 | 0.220 | 0.54 | 0.006 | 79.0
5000 | 0.1 | 0.527 | 1.44 | 0.063 | 95.7
1000 | 0.5 | 0.505 | 0.76 | 0.009 | 60.2
1000 | 0.3 | 0.169 | 0.70 | 0.014 | 94.6
1000 | 0.1 | 0.331 | 3.36 | 0.357 | 98.6
R2SLS (Probit ”Trick”)
5000 | 0.5 | 0.206 | 0.49 | 0.004 | 59.2
5000 | 0.3 | 0.103 | 0.46 | 0.006 | 91.4
5000 | 0.1 | 1.287 | 3.90 | 0.479 | 98.1
1000 | 0.5 | 0.187 | 0.57 | 0.008 | 90.6
1000 | 0.3 | 0.060 | 0.66 | 0.013 | 95.5
1000 | 0.1 | 41.922 | 35.39 | 39.603 | 99.6
Table 5: Simulation results for Complex data setup of Figure 1, based on
b=1000 bootstrapped samples. LPM and Probit indicate first stage models fit
using OLS or a probit regression model respectively. Probit ”Trick” denotes
the method in which first stage predictions of $A_{t}$ using probit are
subsequently used as instruments in 2SLS.

### Full Results: Case Study

Figure 3: Plot of trends of $PP^{Cal}$ and $PP^{t}$ over time. Model |  $\beta_{1}$ (95% CI) |  $\beta_{2}$ (95% CI) |  $\beta_{3}$ (95% CI) |  ATE (95% CI)
---|---|---|---|---
Standard 2SLS $PP^{Cal}$ |  -33.49 (-109.20,40.262) |  35.82 (-70.76,162.19) |  -7.08 (-107.60,84.89) |  -4.75 (-7.11,-2.52)
R2SLS $PP^{Cal}$ |  -9.09 (-51.97,36.11) |  -21.00 (-98.90,34.39) |  -27.61 (-106.90,83.74) |  -2.46 (-4.45,-0.74)
Ridge R2SLS $PP^{Cal}$ |  -1.04 (-1.63,-0.27) |  -0.88 (-1.57,-0.26) |  -0.56 (-1.29,0.07) |  -2.49 (-4.28,-0.65)
Standard 2SLS $PP^{t}$ |  -22.14 (-62.24,8.31) |  0.15 (-55.88,58.43) |  16.64 (-21.18,62.04) |  -5.34 (-8.29,-2.69)
R2SLS $PP^{t}$ |  -18.65 (-44.91,4.91) |  18.39 (-16.18,55.71) |  -2.72 (-30.36,22.62) |  -2.99 (-5.47,-0.81)
Ridge R2SLS $PP^{t}$ |  -1.17 (-1.86,-0.36) |  -0.72 (-1.53,0.01) |  -0.76 (-1.56,-0.02) |  -2.65 (-4.79,-0.66)
Table 6: Full Table of Results for Effect of DPP4 versus SU on HbA1c levels. Characteristics | Correlations
---|---
Time period ($t$) | 0.174*
Age | 0.002
Gender | 0.012
Ethnicity |
White | -0.065
South Asian | 0.094*
Black | -0.043
Other | -0.008
Smoke | -0.082*
HbA1c level (SD) | -0.084*
Median HbA1c | -0.092*
Body Mass Index |
$<25$ | 0.003
$25-29$ | 0.026
$30-34$ | -0.032
$>34$ | 0.022
Years of Diabetes |
$<3$ yrs | 0.029
$3-5$ yrs | -0.036
$>5$ yrs | 0.000
Systolic Blood Pressure |
$<100$ | -0.045
$100-139$ | -0.031
$>140$ | 0.046
Total Serum Cholesterol (mmol/L) (SD) | 0.000
Glomerular Filtration Rate (eGFR) |
$>60$ | 0.016
$30-59$ | -0.025
$<30$ | 0.021
No of Hypoglycemic events | 0.007
Number of strokes | -0.038
Medications (%) |
Statins | 0.025
Beta-Blockers | 0.012
Table 7: Balance Table of correlations between confounders for $PP_{1}^{Cal}$.
Variables with correlations above $0.08$ are marked, indicating potential
imbalance.

## References

  * [1] Hernán M, Robins J. Causal Inference:What If.  Boca Raton: Chapman and Hall/CRC 2020.
  * [2] Baiocchi Michael, Cheng Jing, Small Dylan. Instrumental variable methods for causal inference Statistics in Medicine.  2014;33:4859-4860.
  * [3] Tian Haodong, Patel Ashish, Burgess Stephen. Estimating time-varying exposure effects through continuous-time modelling in Mendelian randomization 2024.
  * [4] Shi Joy, Swanson Sonja, Kraft Peter, Rosner Bernard, Vivo Immaculata, Hernan Miguel. Mendelian Randomization With Repeated Measures of a Time-varying Exposure: An Application of Structural Mean Models Epidemiology.  2022;33:84-94.
  * [5] Cui Y, Michael H, Tanser Frank, Tchetgen E. Instrumental variable estimation of the marginal structural Cox model for time-varying treatments Biometrika.  2021;110.
  * [6] Michael Haben, Cui Yifan, Lorch Scott, Tchetgen Eric. Instrumental Variable Estimation of Marginal Structural Mean Models for Time-Varying Treatment Journal of the American Statistical Association.  2023:1-23.
  * [7] Chen Shuxiao, Zhang Bo. Estimating and Improving Dynamic Treatment Regimes With a Time-Varying Instrumental Variable 2021.
  * [8] Tompsett Daniel, Vansteelandt Stijn, Grieve Richard, Gomes Manuel. Instrumental Variable Approaches For Estimating Time-Varying Treatments in Comparative Effectiveness Research 2024.
  * [9] Robins JM. Marginal structural models versus structural nested models as tools for causal inference. In Statistical models in epidemiology, the environment, and clinical trials.  New York: Springer-Verlag 2000.
  * [10] Angrist Joshua, Imbens Guido. Average Causal Response with Variable Treatment Intensity Journal of the American Statistical Association.  1995;90:431–432.
  * [11] Gudemann Laura, Shields Beverley, Dennis John, Bowden Jack. Just what the doctor ordered: An evaluation of provider preference-based Instrumental Variable methods in observational studies, with application for comparative effectiveness of type 2 diabetes therapy 2023.
  * [12] Sanderson Eleanor, Smith George, Windmeijer Frank, Bowden Jack. An examination of multivariable Mendelian randomization in the single-sample and two-sample summary data settings International Journal of Epidemiology.  2020;49.
  * [13] Burgess Stephen, Foley Christopher, Allara Elias, Staley James, Howson Joanna. A robust and efficient method for Mendelian randomization with hundreds of genetic variants Nature Communications.  2020;11:376.
  * [14] Tian Haodong, Burgess Stephen. Estimation of time-varying causal effects with multivariable Mendelian randomization some cautionary notes International journal of epidemiology.  2023;52.
  * [15] Sanderson Eleanor, Richardson Tom, Morris Tim, Tilling Kate, Smith George. Estimation of causal effects of a time-varying exposure at multiple time points through multivariable mendelian randomization PLOS Genetics.  2022;18:e1010290.
  * [16] Richardson Tom, Sanderson Eleanor, Elsworth Ben, Tilling Kate, Smith George. Use of genetic variation to separate the effects of early and later life adiposity on disease risk: Mendelian randomisation study BMJ.  2020;369.
  * [17] Vansteelandt Stijn, Didelez Vanessa. Robustness and efficiency of covariate adjusted linear instrumental variable estimators Scandinavian Journal of Statistics.  2015;45.
  * [18] Diazordaz Karla, Daniel Rhian, Kreif Noemi. Data-adaptive doubly robust instrumental variable methods for treatment effect heterogeneity Journal de la Societe Francaise de Statistique.  2018\.
  * [19] Diabetes.org.uk https://www.diabetes.org.uk/about-us/about-the-charity/our-strategy/statistics.  Accessed: 2024-29-05.
  * [20] Curtis Helen, Dennis John, Shields Beverley, et al. Time Trends and Geographical Variation in Prescribing of Drugs for Diabetes in England 1998-2017 Diabetes, Obesity and Metabolism.  2018;20.
  * [21] Eriksson Jan, Bodegard Johan, Nathanson David, Thuresson Marcus, Nystrom Thomas, Norhammar Anna. Sulphonylurea compared to DPP-4 inhibitors in combination with metformin carries increased risk of severe hypoglycemia, cardiovascular events, and all-cause mortality Diabetes Research and Clinical Practice.  2016;117.
  * [22] Ling Juan, Cheng Peng, Ge Long, et al. The efficacy and safety of dipeptidyl peptidase-4 inhibitors for type 2 diabetes: a Bayesian network meta-analysis of 58 randomized controlled trials Acta Diabetologica.  2018;56:1-24.
  * [23] Mogstad Magne, Torgovitsky Alexander, Walters Christopher. The Causal Interpretation of Two-Stage Least Squares with Multiple Instrumental Variables American Economic Review.  2021;111:3663-3698.
  * [24] Hernan Miguel, Robins James. Instruments for Causal Inference: An Epidemiologist’s Dream? Epidemiology (Cambridge, Mass.).  2006;17:360-72.
  * [25] Okui Ryo, Small Dylan, Tan Zhiqiang, Robins James. Doubly Robust Instrumental Variable Regression Statistica Sinica.  2012;22.
  * [26] Morris Tim P., White Ian R., Crowther Michael J.. Using simulation studies to evaluate statistical methods Statistics in Medicine.  2019;38:2074-2102.
  * [27] Cragg John G., Donald Stephen. Testing Identifiability and Specification in Instrumental Variable Models Econometric Theory.  1993;9:222-240.
  * [28] Stock James, Yogo Motohiro. Testing for Weak Instruments in Linear IV Regression SSRN eLibrary.  2002;11.
  * [29] Staiger Douglas, Stock James. Instrumental Variables Regression With Weak Instruments Econometrica.  1997;65:557-586.
  * [30] Moler Zapata Silvia, Grieve Richard, Basu Anirban, O’Neill Stephen. How does a local instrumental variable method perform across settings with instruments of differing strengths? A simulation study and an evaluation of emergency surgery Health Economics.  2023;32.
  * [31] Sanderson Eleanor, Windmeijer Frank. A weak instrument F -test in linear IV models with multiple endogenous variables Journal of Econometrics.  2015;190.
  * [32] Bidulka Patrick, Lugo-Palacios David, Carroll Orlagh, et al. Comparative effectiveness of second line oral antidiabetic treatments among people with type 2 diabetes mellitus: emulation of a target trial using routinely collected health data BMJ.  2024;385.
  * [33] Ferrannini Ele, Fonseca V, Zinman B, et al. Fifty-two-week efficacy and safety of vildagliptin vs. Glimepiride in patients with type 2 diabetes mellitus inadequately controlled on metformin monotherapy Diabetes, obesity and metabolism.  2009;11:157-66.
  * [34] Nauck Michael, Meininger Gary, Sheng Dehong, Terranella L, Stein P.P.. Efficacy and safety of the dipeptidyl peptidase-4 inhibitor, sitagliptin, compared with the sulfonylurea, glipizide, in patients with type 2 diabetes inadequately controlled on metformin alone: a randomized, double-blind, non-inferiority trial Diabetes, obesity and metabolism.  2007;9:194-205.
  * [35] Fadini Gian Paolo, Bottigliengo Daniele, D’Angelo Federica, et al. Comparative Effectiveness of DPP-4 Inhibitors Versus Sulfonylurea for the Treatment of Type 2 Diabetes in Routine Clinical Practice: A Retrospective Multicenter Real-World Study Diabetes Therapy.  2018;9.
  * [36] Nagar Shashwat, Napoles Anna, Jordan King, Marino-Ramirez Leonardo. Socioeconomic deprivation and genetic ancestry interact to modify type 2 diabetes ethnic disparities in the United Kingdom EClinicalMedicine.  2021;37:100960.
  * [37] Baltagi Badi H, Griffin James M. A Generalized Error Component Model with Heteroscedastic Disturbances International Economic Review.  1988;29:745-753.
  * [38] Hoerl Arthur, Kennard Robert. Ridge Regression: Biased Estimation for Nonorthogonal Problems Technometrics.  2012;12:55-67.
  * [39] Marquardt Donald W., Snee Ronald D.. Ridge Regression in Practice The American Statistician.  1975;29:3–20.
  * [40] Grant Andrew, Burgess Stephen. Pleiotropy robust methods for multivariable Mendelian randomization 2020.
  * [41] Basu Anirban, Coe Norma, Chapman Cole. 2SLS Versus 2SRI: Appropriate Methods for Rare Outcomes and/or Rare Exposures Health Economics.  2018;27.
  * [42] Angrist Joshua, Pischke Jorn-Steffen. Mostly Harmless Econometrics: An Empiricist’s Companion .  2009\.


* [288] D.C. Moore and A.A. Geraci, _Searching for new physics using optically levitated sensors_ , _Quantum Science and Technology_ 6 (2021) 014008.
  * [289] G. Ranjit, M. Cunningham, K. Casey and A.A. Geraci, _Zeptonewton force sensing with nanospheres in an optical lattice_ , _Phys. Rev. A_ 93 (2016) 053801 [1603.02122].
  * [290] F. Monteiro, S. Ghosh, A.G. Fine and D.C. Moore, _Optical levitation of 10-ng spheres with nano- $g$ acceleration sensitivity_, _Phys. Rev. A_ 96 (2017) 063841.
  * [291] E. Hebestreit, M. Frimmer, R. Reimann and L. Novotny, _Sensing static forces with free-falling nanoparticles_ , _Phys. Rev. Lett._ 121 (2018) 063602.
  * [292] T.M. Hoang, Y. Ma, J. Ahn, J. Bang, F. Robicheaux, Z.-Q. Yin et al., _Torsional optomechanics of a levitated nonspherical nanoparticle_ , _Phys. Rev. Lett._ 117 (2016) 123604.
  * [293] J. Ahn, Z. Xu, J. Bang, Y.-H. Deng, T.M. Hoang, Q. Han et al., _Optically levitated nanodumbbell torsion balance and ghz nanomechanical rotor_ , _Phys. Rev. Lett._ 121 (2018) 033603.
  * [294] R. Reimann, M. Doderer, E. Hebestreit, R. Diehl, M. Frimmer, D. Windey et al., _Ghz rotation of an optically trapped nanoparticle in vacuum_ , _Phys. Rev. Lett._ 121 (2018) 033602.
  * [295] F. Monteiro, S. Ghosh, E.C. van Assendelft and D.C. Moore, _Optical rotation of levitated spheres in high vacuum_ , _Phys. Rev. A_ 97 (2018) 051802.
  * [296] D.C. Moore, A.D. Rider and G. Gratta, _Search for millicharged particles using optically levitated microspheres_ , _Phys. Rev. Lett._ 113 (2014) 251801 [1408.4396].
  * [297] A. Arvanitaki and A.A. Geraci, _Detecting high-frequency gravitational waves with optically levitated sensors_ , _Phys. Rev. Lett._ 110 (2013) 071105.
  * [298] F. Monteiro, G. Afek, D. Carney, G. Krnjaic, J. Wang and D.C. Moore, _Search for composite dark matter with optically levitated sensors_ , 7, 2020.
  * [299] A.A. Geraci, S.B. Papp and J. Kitching, _Short-range force detection using optically cooled levitated microspheres_ , _Phys. Rev. Lett._ 105 (2010) 101101.
  * [300] C.P. Blakemore, A. Fieguth, A. Kawasaki, N. Priel, D. Martin, A.D. Rider et al., _Search for non-newtonian interactions at micrometer scale with a levitated test mass_ , _Phys. Rev. D_ 104 (2021) L061101.
  * [301] D. Carney et al., _Snowmass 2021: Quantum Sensors for HEP Science – Interferometers, Mechanics, Traps, and Clocks_ , in _2022 Snowmass Summer Study_ , 3, 2022 [2203.07250].


# Unsupervised Pre-Training Using Masked Autoencoders for ECG Analysis

Guoxin Wang, Qingyuan Wang, Ganesh Neelakanta Iyer,
Avishek Nag, and Deepu John  This work is partly supported by 1) China
Scholarship Council, 2) JEDAI Project under the Horizon 2020 FET Chist-Era
Program 3) Science Foundation Ireland under Grant number 18/CRT/6183 and 4)
Microelectronic Circuit Centre Ireland. Guoxin Wang, Qingyuan Wang, Avishek
Nag, and Deepu John are with the School of Electrical and Electronic
Engineering, University College Dublin, Dublin 4, Ireland (e-mail:
{guoxin.wang<EMAIL_ADDRESS>{avishek.nag, deepu.john}@ucd.ie;)
and Ganesh Neelakanta Iyer is with Department of Computer Science, National
University of Singapore (e-mail: gni@nus.edu.sg).

###### Abstract

Unsupervised learning methods have become increasingly important in deep
learning due to their demonstrated large utilization of datasets and higher
accuracy in computer vision and natural language processing tasks. There is a
growing trend to extend unsupervised learning methods to other domains, which
helps to utilize a large amount of unlabelled data. This paper proposes an
unsupervised pre-training technique based on masked autoencoder (MAE) for
electrocardiogram (ECG) signals. In addition, we propose a task-specific fine-
tuning to form a complete framework for ECG analysis. The framework is high-
level, universal, and not individually adapted to specific model architectures
or tasks. Experiments are conducted using various model architectures and
large-scale datasets, resulting in an accuracy of 94.39% on the MITDB dataset
for ECG arrhythmia classification task. The result shows a better performance
for the classification of previously unseen data for the proposed approach
compared to fully supervised methods.

###### Index Terms:

Masked Autoencoder, Unsupervised Learning, Big Data, Electrocardiogram

##  I Introduction

Electrocardiogram (ECG) analysis is crucial in diagnosing heart disease and
related biomedical applications [1]. Typically, characteristic features of ECG
such as time intervals, amplitude, and statistical parameters are extracted,
and traditional machine learning methods are used to analyze ECG based on
these features [2]. More recently, deep learning methods have demonstrated
improved performance and effectiveness in biomedical signal analysis [3, 4].
These methods apply supervised learning, and results are often improved by
designing better model architectures. One of the benefits of deep learning is
that it facilitates the extraction of high-dimensional features from the
signal without the need for complex manual pre-processing. In [5], Takalo-
Mattila et al. built an automatic ECG classification system using
Convolutional-Neural-Network (CNN)-based feature extraction. A multilayer
perceptron (MLP) is used to classify ECG beats. This framework achieves an
accuracy of $89.9\%$ when tested with $49712$ samples. Li et al. [6] presented
an arrhythmia classification method that extracted ECG features by Residual
Neural Network (ResNet) and enhanced it by overlapping segmentation method.
They reported an accuracy of $88.9\%$ over $7942$ subjects.

However, traditional supervised learning methods heavily rely on annotated
labels from a single dataset; model training tends to overfit because the
dataset is too small and has possible label errors, limiting the resulting
models’ generalization capability. In addition, many methods claiming
excellent results are performed on intra-patient tasks, a division that
introduces the same record into the training and testing sets, and, therefore,
the high performance that cannot be obtained in real scenarios. Moreover,
these methods often employ engineering techniques that may introduce data
leakage or biases, raising concerns about the credibility of the reported
results. In contrast, unsupervised learning methods offer a compelling
alternative as they do not require labelled data, enabling the utilization of
larger datasets while reducing errors associated with manual annotation.
Masked autoencoder (MAE) is an effective, simple, unsupervised representation
learning strategy proven in computer vision tasks [7], which could be extended
to other research areas [8, 9].

This paper introduces a novel unsupervised pre-training technique based on the
MAE for ECG-related applications and leveraging data augmentation techniques
to improve performance. A task-specific fine-tuning is proposed for downstream
applications. The complete framework is presented systematically, encompassing
all essential aspects ranging from training to testing. To assess its
efficacy, we choose cardiac arrhythmia classification as a case task, and the
framework’s performance is thoroughly evaluated through simulations, providing
valuable information about its capabilities and potential clinical utility.

Figure 1: Overview of the MAE-based unsupervised pre-training technique with
downstream fine-tuning

The contributions of this research are as follows:

  * •

Unsupervised Pre-training and Task-specific Fine-tuning for ECG: This study
presents a novel framework based on MAE for ECG signal analysis. The proposed
framework achieves an accuracy of 94.39% on classifying cardiac arrhythmias in
previously unseen data. Using unsupervised learning techniques, the framework
overcomes the limitations of traditional supervised methods, which require
extensive manual labelling of ECG records.

  * •

Using Larger-Scale datasets: Unlike conventional approaches that rely on
labour-intensive labelling of individual ECG records, the proposed pre-
training reduces the need for independent annotations. This feature
facilitates more accessible and more efficient model training, enabling the
utilization of large-scale datasets without the requirement of extensive
manual labelling. This significantly contributes to improved scalability and
potential for real-world implementations.

The remainder of this paper is organized as follows: Section II details the
proposed unsupervised learning-based technique and fine-tuning. The
experimental results are presented in Section III, and Section IV provides
conclusions and directions for future work.

##  II Method

The complete framework comprises two main parts: pre-training and fine-tuning.
In the pre-training phase, we develop a training strategy for
electrocardiogram signals based on the MAE. We then train the base model using
a sizeable unlabelled dataset. In the fine-tuning phase, we freeze the base
model and train the classifier using a small labelled dataset for the specific
task. An overview of the framework is presented in Figure 1.

###  II-A Pre-train

The predominant approach to classifying cardiac arrhythmias involves
supervised training using a limited dataset. However, this method has a
potential drawback in that the trained model can become over-fitted, resulting
in satisfactory performance on the training dataset but poor generalizability
to other datasets [10]. To alleviate this issue, we propose using an
unsupervised learning method. To this end, MAE-based training has been
devised, described in greater detail below.

The MAE-based solution is conceptually simple in that it removes a portion of
the data and learns to predict what was removed. Its effectiveness has also
been proven in computer vision and natural language processing. Specifically,
this encoder-decoder structure operates by dividing input data into patches,
with the encoder only processing a visible subset of patches, as depicted in
Figure 2b. Subsequently, the decoder reconstructs the input with incomplete
information. Although the reconstructed output may not be perfect, this
approach helps the model better comprehend the input. Once trained, the
decoder can be removed, and the encoder can serve as a practical feature
extractor in other related tasks. An instance of the MAE applied to ECG
signals is shown in Figure 2. MAE is beneficial for using large unlabelled
datasets and can be advantageous for various downstream tasks, particularly
for ECG classification with limited annotated data. In this paper, we utilize
the one-dimensional version of the MAE.

(a) Original

(b) Masked

(c) Reconstructed

Figure 2: An example of MAE, which attempts to reconstruct the original signal
with limited information from the masked signal.

This paper utilizes ConvNeXtV2 [11] to implement a fully convolutional MAE
(FCMAE), which is state-of-the-art in the image classification task. This
method employs a non-symmetric encoder-decoder design and sparse convolution
to reduce computational burden during the pre-training phase. The original
ConvNeXtV2 model was designed for images, whereas our implementation
ConvNeXtV2-1D includes the necessary modifications to accommodate the 1D ECG
signal. The architecture is illustrated in Figure 3.

Figure 3: ConvNeXtV2-1D for ECG applications

###  II-B Fine-tune

####  II-B1 Data Augmentation

Data augmentation enhances the model’s generalization performance in fine-
tuning. Various augmented methods include mixup as described in [12] and
additional white noise. By combining these methods, we ensured that the
fundamental characteristics of the original data, such as the relative
positions of fiducial points and peak intervals, remained unchanged.

####  II-B2 Modelling

We extended the pre-trained model into our framework by adding an MLP head as
a classifier for detecting arrhythmia. We used well-established forward and
backward propagation techniques for supervised learning during this phase.
However, it is essential to note that the pre-trained model remained frozen
throughout this process, ensuring that the previously learned features were
retained without alteration. Therefore, only the newly added classifier
underwent training, allowing it to specialize in accurately identifying
arrhythmia patterns.

##  III Evaluation and Results

###  III-A Datasets

This research used multiple datasets, each serving specific purposes. The
first dataset used for pre-training the model was the PhysioNet / Computing in
Cardiology Challenge 2021 (CINC2021) introduced by Alday et al. [13]. This
large-scale database assembled nine databases with $131,149$ unlabelled
12-lead ECG records for over 1000 hours.

For the subsequent stages of fine-tuning and testing, we employed the MIT-BIH
Arrhythmia Database (MITDB) curated by Moody and Mark [14]. This specific
database consists of 48 records derived from 47 individual subjects. We also
used the St Petersburg INCART 12-lead Arrhythmia Database (INCARTDB) [15] for
fine-tuning. The INCARTDB database includes 75 records from 32 subjects.

To partition the cardiac rhythm data extracted from the MIT-BIH Arrhythmia
Database (MITDB), we adopt a methodology similar to the approach employed by
[16]. The dataset denoted as DS1 is used for fine-tuning, while DS2 serves as
the designated testing dataset. The heartbeats of DS1 and DS2 come from
different individuals. Such division protocol is called in the literature
inter-patient paradigm [17]. Furthermore, DS1 is divided into a ratio of 9:1,
where 90% of the data are allocated for training purposes, while the remaining
10% are allocated for validation. For INCARTDB, we use the whole dataset for
fine-tuning.

The MITDB and INCARTDB database contains labelled annotations for four
classes: N-type (normal beat), SVEB-type (atrial premature beat), VEB-type
(premature ventricular contraction), F-type (fusion of ventricular and normal
beat), and Q-type (unknown beat). Table I summarises the data distribution for
these classes.

TABLE I: data distribution for labelled data

| N | SVEB | VEB | F | Q
---|---|---|---|---|---
INCARTDB | 153563 | 1958 | 20000 | 219 | 6
MITDB-DS1 (Training) | 45298 | 908 | 3597 | 393 | 9
MITDB-DS1 (Validation) | 2392 | 38 | 190 | 21 | 1
MITDB-DS2 | 44225 | 1837 | 3219 | 388 | 7

In the pre-processing step of ECG signals, segmenting them into shorter pieces
using fiducial points is standard practice. Initially, we re-sampled the ECG
signals at a frequency of 360 Hz. Then, we detect all the ’R’ peaks in the ECG
signals for the unlabelled dataset. For labelled datasets, we use annotations,
including peak position and diagnosis information. Next, we extract a segment
of 480 sample points for each ’R’ peak. This segment encompasses 360 points to
the left and 120 points to the right of the ’R’ peak. To ensure consistency
across instances, we normalize each segment to a range between 0 and 1 in the
final.

TABLE IV: Comparison of the proposed system with other published works

Work | Feature Exactor | Classifier | Training Set | Testing Set | Accuracy
---|---|---|---|---|---
Takalo-Mattila et al. [5] | CNN | MLP | MITDB-DS1 | MITDB-DS2 | 89.9%
Li et al. [6] | ResNet | MLP | MITDB-DS1 | MITDB-DS2 | 88.9%
Sellami et al. [18] | CNN | MLP | MITDB-DS1 | MITDB-DS2 | 88.3%
Lin et al. [19] | Morphological Features | Linear Discriminant | MITDB-DS1 | MITDB-DS2 | 91.6%
Asl et al. [20] | Generalized Discriminant Analysis | SVM | MITDB-DS1 + MITDB-DS2 | MITDB-DS1 + MITDB-DS2 | 100%
Chen et al. [21] | Fuducial Features | SVM + MLP | MITDB-DS1 + MITDB-DS2 | MITDB-DS1 + MITDB-DS2 | 100%
Proposed Method | Unsupervised Pre-training | MLP | MITDB-DS1 + INCARTDB | MITDB-DS2 | 94.39%

###  III-B Architecture

Our proposed framework is trained at a high level, and we evaluate it using
different architectures tailored explicitly to the model. To leverage the
simplicity of implementation, we utilize multiple variations of ConvNeXtV2-1D
in the MAE-based training and the supervised training baseline for comparison
and benchmarking purposes. We follow the same configurations of the stage,
block ($B$), and channel ($C$) settings in [11].

  * •

ConvNeXtV2-1D-Atto: $C=40$, $B=(2,2,6,2)$

  * •

ConvNeXtV2-1D-Tiny: $C=96$, $B=(3,3,9,3)$

  * •

ConvNeXtV2-1D-Base: $C=192$, $B=(3,3,27,3)$

###  III-C Experiment Setup

The pre-trained model uses Stochastic Gradient Descent (SGD) with a batch size
of 512 for 500 epochs. Adam optimizer is employed with a batch size of 1024
for 100 epochs for fine-tuning, initializing the learning rate to 0.0003. The
learning rate is gradually reduced using cosine annealing.

To establish a benchmark, we conduct complete supervised training using the
same parameters as the fine-tuning process.

###  III-D Evaluation

MITDB-DS2 is the test set to calculate global performance for different
methods. Accuracy is the main critical metric for the classification task.
Table II shows detailed results of different model architectures and training
strategies. The results show that the proposed method achieved a higher
accuracy than traditional supervised training, which is 90.83% for
ConvNeXtV2-1D-Atto, 91.30% for ConvNeXtV2-1D-Tiny and 91.34% for
ConvNeXtV2-1D-Base.

TABLE II: accuracy with different model architectures and training strategies

| ConvNeXtV2-1D-Atto | ConvNeXtV2-1D-Tiny | ConvNeXtV2-1D-Base
---|---|---|---
Proposed Method | 90.83% | 91.30% | 91.34%
Supervised | 89.48% | 90.58% | 88.54%

In addition, Table III shows that adding different datasets for fine-tuning
and supervised learning is helpful for performance improvement. MAE-based
training with MITDB-DS1 and INCARTDB fine-tuning achieved higher accuracy
among multiple architecture complexity, which is 94.39% for
ConvNeXtV2-1D-Atto, 93.98% for ConvNeXtV2-1D-Tiny and 93.89% for
ConvNeXtV2-1D-Base.

TABLE III: accuracy with different fine-tuning datasets

| ConvNeXtV2-1D-Atto | ConvNeXtV2-1D-Tiny | ConvNeXtV2-1D-Base
---|---|---|---
Proposed Method (MITDB-DS1) | 90.83% | 91.30% | 91.34%
Proposed Method (MITDB-DS1, INCARTDB) | 94.39% | 93.98% | 93.89%

Table IV compares the proposed framework and the current, reliable state of
the art. [5, 6, 18], who have reported performance on the MITDB-DS2 with
supervised deep learning methods, which is less than ours. Furthermore, their
system used a single dataset for training and did not take full advantage of
the available ECG dataset, resulting in low accuracy. Furthermore, [19] used
methods based on machine learning and required a lot of manual handling but
also reported low precision. [20, 21] have reported an accuracy of 100%, but
they actually conducted the intra-patient test, where heartbeats of the same
records probably appear in training and the testing dataset. However, in a
realistic scenario, a fully automatic method will find patients’ heartbeats
different from those they used to learn in the training phase; hence, the high
accuracy they report is questionable. Compared with these frameworks, our
proposed framework achieves higher accuracy within the same task, fully uses
existing datasets, and meets practical needs.

## IV Conclusion

Our study introduced a new MAE-based cardiac arrhythmia classification system.
The system uses unsupervised learning to learn generic ECG information and
classify arrhythmia after fine-tuning. Experiments show that the proposed
approach improves performance compared to traditional methods. Future work
includes using different unsupervised learning approaches, more architectures,
larger datasets, model compression, embedded system deployment, and transfer
learning exploration.

## References

  * [1] Eduardo José da S. Luz, William Robson Schwartz, Guillermo Cámara-Chávez, and David Menotti.  ECG-based heartbeat classification for arrhythmia detection: A survey.  Comput. Methods Programs Biomed., 127:144–164, April 2016.
  * [2] Shenda Hong et al.  ENCASE: an ENsemble ClASsifiEr for ECG Classification Using Expert Features and Deep Neural Networks.  ResearchGate, September 2017.
  * [3] Li Xiaolin, Hans Vandierendonck, Dimitrios S. Nikolopoulos, Bo Ji, Barry Cardiff, and Deepu John.  Decentralised biomedical signal classification using early exits.  In 2023 21st IEEE Interregional NEWCAS Conference (NEWCAS), pages 1–2, 2023.
  * [4] Gawsalyan Sivapalan et al.  Interpretable rule mining for real-time ecg anomaly detection in iot edge sensors.  IEEE Internet of Things Journal, 10(15):13095–13108, 2023.
  * [5] Janne Takalo-Mattila et al.  Inter-Patient ECG Classification Using Deep Convolutional Neural Networks.  IEEE Computer Society, August 2018.
  * [6] Yuanlu Li, Renfei Qian, and Kun Li.  Inter-patient arrhythmia classification with improved deep residual convolutional neural network.  Comput. Methods Programs Biomed., 214:106582, February 2022.
  * [7] Kaiming He et al.  Masked autoencoders are scalable vision learners.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 16000–16009, 2022.
  * [8] Guoxin Wang et al.  Low complexity ecg biometric authentication for iot edge devices.  In 2020 IEEE International Conference on Integrated Circuits, Technologies and Applications (ICTA), pages 145–146, 2020.
  * [9] Conor Smyth et al.  Continuous user authentication using iot wearable sensors.  In 2021 IEEE International Symposium on Circuits and Systems (ISCAS), pages 1–5, 2021.
  * [10] Xue Ying.  An Overview of Overfitting and its Solutions.  J. Phys. Conf. Ser., 1168(2):022022, February 2019.
  * [11] Sanghyun Woo, Shoubhik Debnath, Ronghang Hu, Xinlei Chen, Zhuang Liu, In So Kweon, and Saining Xie.  Convnext v2: Co-designing and scaling convnets with masked autoencoders, 2023.
  * [12] Hongyi Zhang, Moustapha Cisse, Yann N. Dauphin, and David Lopez-Paz.  mixup: Beyond Empirical Risk Minimization.  arXiv, October 2017.
  * [13] Erick A. Perez Alday et al.  Classification of 12-lead ECGs: the PhysioNet/Computing in Cardiology Challenge 2020.  Physiol. Meas., 41(12):124003, December 2020.
  * [14] George B. Moody and Roger G. Mark.  Mit-bih arrhythmia database.  https://doi.org/10.13026/C2F305, 1992.
  * [15] V. Tihonenko, A. Khaustov, S. Ivanov, and A. Rivin.  St.-petersburg institute of cardiological technics 12-lead arrhythmia database.  https://doi.org/10.13026/C2V88N, 2007.
  * [16] Philip de Chazal, M. O’Dwyer, and R.B. Reilly.  Automatic classification of heartbeats using ecg morphology and heartbeat interval features.  IEEE Transactions on Biomedical Engineering, 51(7):1196–1206, 2004\.
  * [17] Gaël de Lannoy, Damien Francois, Jean Delbeke, and Michel Verleysen.  Weighted conditional random fields for supervised interpatient heartbeat classification.  IEEE Transactions on Biomedical Engineering, 59(1):241–247, 2012\.
  * [18] Ali Sellami and Heasoo Hwang.  A robust deep convolutional neural network with batch-weighted loss for heartbeat classification.  Expert Syst. Appl., 122:75–84, May 2019.
  * [19] Chun-Cheng Lin and Chun-Min Yang.  Heartbeat Classification Using Normalized RR Intervals and Morphological Features.  Math. Prob. Eng., 2014, May 2014.
  * [20] Babak Mohammadzadeh Asl, Seyed Kamaledin Setarehdan, and Maryam Mohebbi.  Support vector machine-based arrhythmia classification using reduced features of heart rate variability signal.  Artif. Intell. Med., 44(1):51–64, September 2008.
  * [21] Huan Chen, Bo-Chao Cheng, Guo-Tan Liao, and Ting-Chun Kuo.  Hybrid classification engine for cardiac arrhythmia cloud service in elderly healthcare management.  Journal of Visual Languages & Computing, 25(6):745–753, December 2014.