mainReferences appendixAppendix References

# Estimating time-varying exposure effects through continuous-time modelling
in Mendelian randomization

Haodong Tian1, Ashish Patel1, Stephen Burgess1,2
1 MRC Biostatistics Unit, School of Clinical Medicine,
University of Cambridge, Cambridge, UK
2 British Heart Foundation Cardiovascular Epidemiology Unit,
Department of Public Health and Primary Care,
University of Cambridge, Cambridge, UK haodong.tian@mrc-
bsu.cam.ac.ukashish.patel<EMAIL_ADDRESS>
###### Abstract

Mendelian randomization is an instrumental variable method that utilizes
genetic information to investigate the causal effect of a modifiable exposure
on an outcome. In most cases, the exposure changes over time. Understanding
the time-varying causal effect of the exposure can yield detailed insights
into mechanistic effects and the potential impact of public health
interventions. Recently, a growing number of Mendelian randomization studies
have attempted to explore time-varying causal effects. However, the proposed
approaches oversimplify temporal information and rely on overly restrictive
structural assumptions, limiting their reliability in addressing time-varying
causal problems. This paper considers a novel approach to estimate time-
varying effects through continuous-time modelling by combining functional
principal component analysis and weak-instrument-robust techniques. Our method
effectively utilizes available data without making strong structural
assumptions and can be applied in general settings where the exposure
measurements occur at different timepoints for different individuals. We
demonstrate through simulations that our proposed method performs well in
estimating time-varying effects and provides reliable inference results when
the time-varying effect form is correctly specified. The method could
theoretically be used to estimate arbitrarily complex time-varying effects.
However, there is a trade-off between model complexity and instrument
strength. Estimating complex time-varying effects requires instruments that
are unrealistically strong. We illustrate the application of this method in a
case study examining the time-varying effects of systolic blood pressure on
urea levels.

Keywords: instrumental variables, genetics, functional data analysis,
principal components, instrument strength, identification-robust inference

##  1 Introduction

Mendelian randomization (MR) utilizes genetic variants as instrumental
variables (IVs) to assess the causal effects of modifiable exposures on health
outcomes [10, 21]. In reality, most modifiable exposures are subject to change
over time, indicating that the exposure effects on outcomes should be modelled
over time to incorporate time-varying patterns. Examining time-varying effects
can provide valuable insights into underlying mechanisms. For example, such
analyses can provide evidence on the critical time period of treatment that
has greatest effect on the outcome, or the long-term effect of some risk
factors. Time-varying causal studies have recently emerged as a popular
research area, leading to the development of numerous models and methods
tailored to time-related questions. A recent review of time-varying causal
models and methods can be found in the paper [26].

Though time-varying causal models and methods are numerous, they are hard to
tailor to MR analysis. Many of the methods (like those based on structural
nested mean models [31, 32] or marginal structural mean models [33, 22]) rely
on detailed individual-level data (e.g., comprehensive longitudinal
measurements in one sample) and strong assumptions [43], and therefore, they
may not be generally applicable to common real datasets. Most existing MR
approaches for time-varying analysis consider the exposure at specific time
points or overlook the time factor, possibly failing to provide a reliable
understanding of the time-varying effects. Since the genetic variants are
randomized at conception and their effects on the exposure are long-term
rather than short-term fluctuations [40], MR estimates are usually interpreted
as lifelong effects [7, ch. 6.1.1]. However, the lifelong effect usually
represents a cumulative effect and cannot provide time-varying information
during particular time periods. The challenges in time-varying MR studies have
been discussed in a recent paper [27].

In recent years, a growing body of literature has focused on the exploration
of time-varying effects using MR. This includes work that considers repeated
exposure measurements at different time points as separate exposures in a
multivariable MR (MVMR) model [30, 39, 29, 37]. However, these methods rely on
a model represented by a directed acyclic graph (DAG) where each node
corresponds to the phenotype at a specific timepoint—potentially predetermined
by researchers. However, there is an absence of a clear and reasonable data-
generating mechanism supporting the models in the time-varying context. Using
DAGs with discrete time points for time-varying causal inference has been
criticized [2, 1]. Additionally, there is difficulty in interpreting the
effects from such models. In particular, it is not clear what intervention
corresponds to the estimated parameters. A recent paper has illustrated that
employing the same risk factor measured at a limited number of measurement
times in MVMR can potentially lead to misleading results [45].

Another line of research has explored the interpretation of MR estimates with
a time-varying exposure. For instance, recent work suggests that MR with a
single exposure measurement can yield misleading estimates if interpreted as
the effect of varying the exposure in a mature population [18, 19].
Furthermore, in scenarios with multiple measurement timepoints for the
exposure, MR may fail to provide a solid basis for reasonable causal null
testing over certain time regions [42]. These challenges arise from the
limited trajectory information considered for exposures, which violates the IV
exclusion restriction assumption as the exposure continuously exists over
time.

When we have limited exposure information available, such as measurements at
specific time points, MR can still serve as a valid method for estimating
specific cumulative time-varying quantities. For example, when the genetic
effects remain constant, MR results can be interpreted as representing the
cumulative effect over an individual’s lifetime [18]. However, for more
detailed time-varying conclusions, we rely on changes over time in the effects
of the genetic variants on the exposure, which are referred to as time-varying
genetic-exposure associations [7, ch. 8]. Some evidence suggests that genetic
associations with certain exposures can indeed vary over time [20], which
bolsters the feasibility of studying time-varying effects.

Given that time is continuous, and the corresponding effect of an exposure
over time is also continuous, we should model the time-varying effect as a
functional objective and incorporate the time-varying exposure trajectory into
the model. Since the functional objective is infinite-dimensional, we utilize
functional principal component analysis to achieve dimension reduction by
extracting low-dimensional factors from exposure trajectories that can
sufficiently explain time-varying information. Then, we integrate these low-
dimensional factors into a multivariable MR model. These multivariable MR
estimates are used to back out the functional time-varying exposure effect.

The aim of this work is to explore a method for time-varying MR analysis that
can be applied in general settings. Some key features of our method are: we
treat the time factor and the effect defined over time as a functional
objective and use continuous-time modelling to study the time-varying effect;
and we do not rely on strong parametric assumptions for the exposure
generation mechanism. We allow the exposure to be measured at different time
points for every individual, and each individual can have sparse exposure
measurements. We also consider practical issues relating to time-varying MR
studies under weak instruments.

The manuscript is structured as follows. In Section 2, we introduce general
time-varying models and discuss the corresponding IV assumptions. We also
delve into the issues related to identification under different settings and
the challenges surrounding identification when employing existing MR methods.
Section 3 demonstrates the method of using functional principal component
analysis to build a multivariable MR model and estimate time-varying effects
as a functional objective. In Section 4, we present simulation examples with
various time-varying instrument-exposure and exposure-outcome model scenarios.
In Section 5, we provide an illustrative example applying our methodology to
real data to study the effect of blood pressure on urea levels. We conclude in
Section 6 with discussions. The code implementing our time-varying MR analysis
is available at: https://github.com/HDTian/TVMR.

##  2 Models

We denote the time-varying exposure as a stochastic process starting from
timepoint $s$ by $\vec{X}^{s}=\\{X(t);t\geq s\\}$, where $X(t)$ represents the
exposure level at timepoint $t$. We let $\vec{X}=\vec{X}^{0}$. Furthermore,
$Z$ and $U(t)$ represent a genetic variant and (possibly time-varying)
unobserved confounders at timepoint $t$, respectively. The outcome of interest
at timepoint $T$ is denoted as $Y_{T}$. The graphical representation of the
relationships between these variables is illustrated in Figure 1, which can be
viewed as the continuous time version of the traditional instrumental variable
DAG. We denote the potential outcome resulting from a complete exposure
trajectory intervention $\vec{x}$ and an instrument intervention $z$ as
$Y_{T}(\vec{x},z)$. We assume that all potential outcomes are well-defined,
adhering to the Stable Unit Treatment Value Assumption (SUTVA) [35].

To investigate time-varying effects over the interval $[0,T]$ using
instrumental variable (IV) methods, we require that any genetic variant used
as an instrument satisfies the following three core IV assumptions:

  * i.

(Relevance) $cov(Z,X(t))\neq 0$ for some $t$ over $[0,T]$

  * ii.

(Exchangeability) $Z\rotatebox[origin={c}]{90.0}{$\models$}Y_{T}(\vec{x},z)$
for any $\vec{x}$ and $z$

  * iii.

(Exclusion restriction)
$Y_{T}(\vec{x},z_{1})=Y_{T}(\vec{x},z_{2})=Y_{T}(\vec{x})$ for any $z_{1}$,
$z_{2}$ and $\vec{x}$

Those assumptions are based on the potential outcome framework [34] and are
analogous to the classical IV assumptions [16]. However, in the time-varying
case, the single exposure is replaced by the exposure trajectory, which is a
functional variable depending on time. Relevance implies that each instrument
should be associated with the exposure over some time regions. The exclusion
restriction assumption states that each instrument can only affect the outcome
through the exposure. Note that we do not require the exchangeability
additionally satisfying $Z\rotatebox[origin={c}]{90.0}{$\models$}\vec{X}(z)$,
where $\vec{X}(z)$ represents the counterfactual exposure trajectory given the
instrument level $z$. This means that the instrument can serve as a surrogate
rather than a causal IV (e.g., a variant may be in correlation with some valid
causal variants) [14].

The three IV core assumptions are sufficient to test the joint sharp null
hypothesis $H_{0}:\,Y_{T}(\vec{x}_{1})-Y_{T}(\vec{x}_{2})=0$ for any
individuals and $\vec{x}_{1}$ and $\vec{x}_{2}$, which means that the exposure
has no effect on the outcome. However, such a conclusion cannot be further
specified for specific time sub-regions. For example, to test the local sharp
null hypothesis that $Y_{T}(\vec{x}^{s}_{1})-Y_{T}(\vec{x}^{s}_{2})=0$ for an
intermediate time point $0<s<t<T$, the exclusion restriction originally
satisfied for $[0,T]$ may be violated for $[s,T]$ (e.g., there may be a path
$Z\to X_{0}\to Y_{T}$ in Figure 1). When the joint sharp null hypothesis is
rejected, it is difficult to gain further insight on how a treatment effect
may vary over time.

Figure 1: The graph illustrates the IV assumptions for the time-varying
Mendelian randomization setting. $X_{0}$, $X_{t}$, and $X_{T}$ represent the
exposure at the time points $0$, $t$, and $T$, respectively, where $t$ is an
arbitrary intermediate time point. $Z$ represents the time-invariant
instrument existing at time point $0$ (e.g., the genetic variants).
$\vec{U}=\\{U(t);0\leq t\leq T\\}$ represents the process of the confounders.
$Y_{T}$ is the outcome of interest, respectively. There should be no arrows
from the exposure (or confounders) at a later time point to the outcome at
present or past time points, as the future cannot influence the past.

To investigate the time-varying effects beyond testing the joint sharp null,
additional assumptions regarding the model are necessary. We assume the
outcome is defined at the end timepoint $T$, and the structural equation for
the outcome $Y$ is as follows:

$Y=\beta_{0}+\int_{0}^{T}\beta(t)X(t)dt+g_{Y}(\vec{U},\epsilon_{Y})$ (1)

Here, $g_{Y}(\cdot)$ is an unknown function, $\vec{U}=\\{U(t);0\leq t\leq
T\\}$ is an unmeasured confounder trajectory, consisting of time-invariant
variables and stochastic processes from the past up to the end time point. The
term $\epsilon_{Y}$ denotes an exogenous error, and $\beta(t)$ represents the
direct effect of the exposure at time point $t$ on the outcome, corresponding
to the arrow from $X_{t}$ to $Y_{T}$ in Figure 1. The function $\beta(t)$
provides complete information about the time-varying effect of the exposure on
the outcome. Under any known time-varying treatment $\vec{a}$, one can express
the effect of the exposure on the outcome via the counterfactual contrast

$\mathbb{E}(Y(\vec{X}+\vec{a})-Y(\vec{X}))=\int_{0}^{T}\beta(t)a(t)dt$ (2)

We treat $\beta(t)$ as the functional objective of interest and call it the
effect function.

We note that the discrete-time models used in recent studies [30, 39, 29, 37]
can be regarded as a special case of our time-varying model. For example,
assuming $\beta(t)=\sum_{k=1}^{K}\beta_{k}\delta_{D}(t-t_{k})$ for certain
time points ${t_{k},k=1,\ldots,K}$, where $\delta_{D}(t)$ represents the Dirac
delta function.

We now take the inverse-variance-weighting (IVW) method as a simple estimation
method, although other methods like the two-stage least squares (2SLS) lead to
the same conclusion. Given the core IV assumptions, we can fit the following
weighted regression: (see Supplementary Text S1)

$\hat{\theta}_{j}=\int_{0}^{T}\beta(t)\hat{\alpha}_{j}(t)dt+\epsilon_{j}\qquad\epsilon_{j}\sim\mathcal{N}(0,s.e.(\hat{\theta}_{j})^{2})\qquad
j=1,2\ldots$ (3)

where $\hat{\theta}_{j}$ represents the estimated association between the
$j$-th instrument and the outcome, while $\hat{\alpha}_{j}(t)$ represents the
estimated association between the $j$-th instrument and the exposure at time
point $t$. Equation (3) captures the fundamental idea of time-varying MR
studies: the identification of the time-varying effect $\beta(t)$ relies on
the changes in exposure over time induced by the genetic variants (i.e.,
$\alpha_{j}(t)$) [7, ch. 8]. If $\alpha_{j}(t)$ remains constant for any time
point for all instruments, it becomes impossible to study the shape of
$\beta(t)$. For instance, Labrecque et al. demonstrated that in such cases,
only the cumulative effect can be identified [18]. However, identifying the
functional objective $\beta(t)$ necessitates some variation in shape among
$\\{\alpha_{j}(t),j=1,2,\ldots\\}$ (i.e. the time-varying genetic-exposure
association across the genetic instruments). We assume that the
$\alpha_{j}(t)$ are sufficiently variable to allow estimation of the time-
varying effect $\beta(t)$. A sufficient condition is that there is no
timepoint when the $\alpha_{j}(t)$ are all zero, and there is no time period
$t\in[t_{1},t_{2}]$ where the $\alpha_{j}(t)$ are all constant or all linearly
dependent.

An additional challenge arises from the high uncertainty typically associated
with estimating the functions $\hat{\alpha}_{j}(t)$. This issue is
particularly prominent in time-varying studies where data are limited for
certain time periods. To overcome this limitation and enhance power in time-
related analyses, information may be borrowed across time.

##  3 Methods

Assume we can express any individual’s exposure trajectory in a linear
additive form as

$X(t)=\mu(t)+\sum_{k=1}^{K}\xi_{k}\phi_{k}(t)+e(t)$ (4)

where the functional terms $\\{\mu(t),\phi_{k}(t)\\}$ are known time-varying
trends that are common for all individuals; the nonfunctional terms
$\\{\xi_{k}\\}$ are the individual-specific variables, which summarize the
variation in the exposure over time through $\phi_{k}(t)$; $e(t)$ is a time-
varying error term for the exposure and assumed to be sufficiently small so
that the structural model (1) becomes

$\begin{split}Y_{i}=&\beta_{0}+\int_{0}^{T}\beta(t)X_{i}(t)dt+g_{Y}(U_{i},\epsilon_{Y,i})\\\
=&\beta_{0}+\int_{0}^{T}\beta(t)\left[\mu(t)+\sum_{k=1}^{K}\xi_{k,i}\phi_{k}(t)+e_{i}(t)\right]dt+g_{Y}(U_{i},\epsilon_{Y,i})\\\
=&\underbrace{\left[\beta_{0}+\int_{0}^{T}\beta(t)\mu(t)dt\right]}_{\beta_{0}^{\ast}}+\sum_{k=1}^{K}\underbrace{\left[\int_{0}^{T}\beta(t)\phi_{k}(t)dt\right]}_{\beta_{k}^{\ast}}\xi_{k,i}+\underbrace{\int_{0}^{T}\beta(t)e_{i}(t)dt+g_{Y}(U_{i},\epsilon_{Y,i})}_{\epsilon_{i}}\end{split}$
(5)

which is the classical form of MVMR where $\\{\xi_{k}\\}$ are considered the
exposure variables, and their effects are given by $\\{\beta^{\star}_{k}\\}$.
We refer to $\\{\xi_{k}\\}$ as pseudo-exposures.

Using the linear additive form for the exposure trajectory reduces the
original infinite-dimensional problem with respect to the effect function
$\beta(t)$ to a $K$-dimensional problem with respect to the components
$\\{\xi_{k}\\}$ that can be addressed through IV analysis. The remainder term
$\int_{0}^{T}\beta(t)e_{i}(t)dt$ may lead to a violation of the IV core
assumptions, especially when the trajectory approximation is poor. We refer to
this as structural pleiotropy, which can be negligible if $e(t)$ is
sufficiently small or is independent of the instruments.

Several methods exist to accomplish such dimension reduction. One approach
involves making a parametric assumption for individual exposure trajectories
and extrapolating or interpolating the exposure for each individual.
Alternatively, one can first estimate the genetic association with the
exposure $\hat{\alpha}(t)$ at multiple timepoints and then smooth these values
over time. However, these methods come with their own set of challenges. The
former method relies on parametric assumptions. The latter method is likely to
suffer from high estimate uncertainty for any specific timepoint due to the
limited sample size for the exposure at certain measured timepoints. This
motivates us to consider more efficient dimension reduction methods to avoid
these problems.

###  3.1 Estimating $\beta(t)$

Consider the functional approximation (5). In the case that
$cov(\xi_{k},\epsilon)=0$ for any $k$, we can estimate the effect shape
$\beta(t)$ and its confidence intervals by the following two steps:

  * 1.

First, fit linear regression of $Y$ on the pseudo-exposures $\\{\xi_{k}\\}$ to
obtain the estimates $\hat{\boldsymbol{\beta}}^{\ast}$ and the covariance
matrix $\Sigma_{\boldsymbol{\beta}^{\ast}}$.

  * 2.

Invert the equation $\beta_{k}^{\ast}=\int\beta(t)\phi_{k}(t)dt$ using the
functions $\phi(t)$ and plug in the estimates $\hat{\beta}^{\ast}_{k}$ (as
well as $\Sigma_{\hat{\boldsymbol{\beta}}^{\ast}}$) to obtain
$\hat{\beta}(t)$. Section 3.2.3 describes further details of inverting the
integral equation $\beta_{k}^{\ast}=\int\beta(t)\phi_{k}(t)dt$.

Since the exposure trajectory $\\{X(t)\\}$ (and thus $\\{\xi_{k}\\}$) can be
confounded with the outcome $Y$, $cov(\xi_{k},U)\neq 0$ and
$cov(\xi_{k},\epsilon)=0$ may not hold, and therefore we instrument the
pseudo-exposures $\xi_{k}$ using genetic variants.

###  3.2 Multi-principal-components Mendelian randomization

In this section, we will introduce (i) how to apply data-driven methods for
functional dimension reduction, (ii) using instruments to study time-varying
effects, and finally, (iii) how to derive the effect function estimator.

####  3.2.1 Functional principal component analysis

One effective approach to achieve dimension reduction for the exposure
trajectories and reduce the variance of the remainder term $e(t)$ is called
functional principal component analysis (FPCA). The Karhunen–Loève theorem
states that, under some regularity conditions,
$X(t)=\mu(t)+\sum_{k=1}^{\infty}\xi_{k}\phi_{k}(t)$ where
$\mu(t)=\mathbb{E}(X(t))$, $\xi_{k}$ represents the principal component (PC)
value with respect to the $k$-th eigenfunction $\phi_{k}(t)$, and
$\lambda_{k}:=\text{var}(\xi_{k})$ denotes the eigenvalues such that
$\lambda_{1}>\lambda_{2}>\cdots$. FPCA transforms the variable trajectories
into a PC-weighted sum of eigenfunctions, ordered by the primary PC scores
that explain more variance. For more details, please refer to Supplementary
Text S2. Therefore, we can approximate the exposure trajectory using the first
$K$ PCs that explain a sufficient amount of variation, making $e_{i}(t)$
sufficiently small. In practice, the choice of $K$ can be based on the number
of principal components required to explain a certain proportion of exposure
variation. Since the available data may have only sparse measurement time
points and the exposure measurements may be observed at different time points
for different individuals, one method for conducting FPCA in such scenarios is
principal component analysis through conditional expectation (PACE) [47, 46].
PACE helps to estimate individual PC scores $\xi_{i,k}$ and the eigenfunctions
$\phi_{k}(t)$ based on the properties of FPCA. Further details can be found in
Supplementary Text S3.

The analogous diagram illustrating Multi-principal-components Mendelian
randomization (MPCMR) is provided in Figure 2. In this graph, the pseudo-
exposures $\\{\xi_{1},\ldots,\xi_{K}\\}$ can be interpreted as latent time-
invariant variables mediating the path from the genetic variants to the time-
varying exposure [23, 37].

Figure 2: The graph of a Mendelian randomization analysis under the MPCMR
model. $Z$ represents the high-dimensional genetic instruments,
${\xi_{1},\ldots,\xi_{K}}$ are the $K$ principal components associated with
the individual exposure curve ${X(t)}$, $\vec{U}=\\{U(t);0\leq t\leq T\\}$
stands for unmeasured confounders, and $Y$ denotes the outcome.

####  3.2.2 MPCMR fitting

Once we have computed the principal components (PCs), we can fit the MPCMR
models as described in (5) using standard IV methods. We can estimate the
effect function using the same two-step procedure as in the association model
described in Section 3.1. We consider the continuously updating generalized
method of moments (GMM) [13] rather than the IVW method for estimation, as it
is less sensitive to weak instruments [4]. To proceed, assume that all
variables (including the genetic variants, principal components, and the
outcome) have been mean centered. We denote the sample mean of any random
variable $X$ with a sample size of $n$ as $\hat{\mathbb{E}}_{n}(X)$. We define
the $J$-vector of instruments as $\boldsymbol{Z}$ and the $K$-vector of
pseudo-exposures (principal components) as $\boldsymbol{\xi}$. The GMM
estimator is:

$\hat{\boldsymbol{\beta}}^{\ast}=\operatorname*{arg\,min}_{\boldsymbol{\beta}}\hat{\boldsymbol{g}}(\boldsymbol{\beta})\hat{\Omega}(\boldsymbol{\beta})^{-1}\hat{\boldsymbol{g}}(\boldsymbol{\beta})$
(6)

where
$\hat{\boldsymbol{g}}({\boldsymbol{\beta}})=\hat{\mathbb{E}}_{n}(\boldsymbol{Z}(Y-\boldsymbol{\xi}^{T}{\boldsymbol{\beta}}))$
and
$\hat{\Omega}(\boldsymbol{\beta})=\hat{\mathbb{E}}_{n}(\boldsymbol{Z}\boldsymbol{Z}^{T})\hat{\mathbb{E}}_{n}((Y-\boldsymbol{\xi}_{i}^{T}{\boldsymbol{\beta}})^{2})$.
Under the assumption that the covariance between the composite error term
$\epsilon$ in Equation (5) and the instruments $\boldsymbol{Z}$ is assumed to
be zero, homoskedastic errors, and standard GMM regularity conditions [24], we
have

$\sqrt{n}(\hat{\boldsymbol{\beta}}^{\ast}-\boldsymbol{\beta}^{\ast})\overset{D}{\to}\mathcal{N}(\boldsymbol{0},(\boldsymbol{G}^{T}\Omega(\boldsymbol{\beta}^{\ast})^{-1}\boldsymbol{G})^{-1})$
(7)

where $\boldsymbol{G}=\mathbb{E}(-\boldsymbol{Z}\boldsymbol{\xi}^{T})$ and
$\boldsymbol{\Omega}(\boldsymbol{\beta}^{\ast})=\operatorname*{plim}_{n\rightarrow\infty}\hat{\Omega}(\boldsymbol{\beta}^{\ast})$.
We can obtain the estimated variance-covariance matrix of
$\hat{\boldsymbol{\beta}}^{\ast}$, denoted by $\hat{\Sigma}$, with plug-in
estimates $\hat{\boldsymbol{\beta}}^{\ast}$ and
$\hat{\boldsymbol{G}}=\hat{\mathbb{E}}_{n}(-\boldsymbol{Z}\boldsymbol{\xi}^{T})$.
That is,
$\hat{\Sigma}=(\hat{\boldsymbol{G}}^{T}\hat{\Omega}(\hat{\boldsymbol{\beta}}^{\ast})^{-1}\hat{\boldsymbol{G}})^{-1}$.

####  3.2.3 Curve inversion

We can proceed to invert the equation
$\beta_{k}^{\ast}=\int\beta(t)\phi_{k}(t)dt$ to obtain the functional
objective $\beta(t)$. Define
$\boldsymbol{\phi}=(\phi_{1}(t)\,\phi_{2}(t)\,\ldots\,\phi_{K}(t))^{T}$. Given
that $\boldsymbol{\beta}^{\ast}=\int\boldsymbol{\phi}(t)\beta(t)dt$ and the
eigenfunction property
$\int\boldsymbol{\phi}(t)\boldsymbol{\phi}^{T}(t)dt=\boldsymbol{I}$, we can
immediately find one solution for the equation, which is
$\beta(t)=\boldsymbol{\phi}^{T}(t)\mathbb{\boldsymbol{\beta}^{\ast}}$.
Therefore, we have one final MR estimator for the time-varying effect
function:

$\hat{\beta}(t)=\boldsymbol{\phi}^{T}(t)\hat{\boldsymbol{\beta}}^{\ast}=\sum_{k=1}^{K}\phi_{k}(t)\hat{\beta}^{\ast}_{k}$
(8)

and the pointwise standard errors
$\hat{\sigma}(t)=\sqrt{\boldsymbol{\phi}^{T}(t)\hat{\Sigma}\boldsymbol{\phi}(t)}$.

Note that the inversion of $\boldsymbol{\beta}^{\ast}$ presents an ill-posed
inverse problem. This means that the time-varying effect $\beta(t)$ may not be
identified even if we have a consistent estimator of
${\boldsymbol{\beta}}^{\ast}$ [15]. There exist infinitely many solutions for
$\beta(t)$ such that the equation
$\boldsymbol{\beta}^{\ast}=\int\boldsymbol{\phi}(t)\beta(t)dt$ hold. To
address this issue, additional parametric assumptions may be useful to study
the objective function [25]. These assumptions can help constrain the possible
solutions and improve the identifiability of the time-varying effect.

In the general case, we assume that the time-varying function can be
represented as a linear additive combination of basis functions, such that
$\beta(t)=\sum_{l=1}^{L}\gamma_{l}b_{l}(t)$. Here, $\\{b_{l}(t)\\}$ are known
basis functions, and $\\{\gamma_{l}\\}$ are the parameters of interest
representing the contribution of each basis function to the time-varying
effect. The structural outcome equation can be rewritten to accommodate this
parametric form as follows:

$\begin{split}Y_{i}=&\beta_{0}+\sum_{k=1}^{K}\left[\int_{0}^{T}\beta(t)\phi_{k}(t)dt\right]\xi_{k,i}+\epsilon_{i}\\\
=&\beta_{0}+\sum_{k=1}^{K}\left[\int_{0}^{T}\left[\sum_{l=1}^{L}\gamma_{l}b_{l}(t)\right]\phi_{k}(t)dt\right]\xi_{k,i}+\epsilon_{i}\\\
=&\beta_{0}+\sum_{l=1}^{L}\gamma_{l}\underbrace{\left[\sum_{k=1}^{K}\int_{0}^{T}b_{l}(t)\phi_{k}(t)dt\,\xi_{k,i}\right]}_{\xi^{\ast}_{l,i}}+\epsilon_{i}\end{split}$
(9)

where we refer to $\boldsymbol{\xi}^{\ast}\in\mathbb{R}^{L\times 1}$ as the
transformed exposure (principal components) vector, which satisfies
$\boldsymbol{\xi}^{\ast}=\boldsymbol{B}^{T}\boldsymbol{\xi}$, where
$\boldsymbol{B}=\int_{0}^{T}\boldsymbol{\phi}(t)\boldsymbol{b}^{T}(t)dt\in\mathbb{R}^{K\times
L}$ is a known transforming matrix. Therefore, we can fit the MPCMR with the
exact same procedures, but replace the exposures with the transformed
exposures. In GMM, we can simply replace the $K$-vector of pseudo-exposure
with the $L$-vector of transformed exposure. Note that $L\leq K$ is required
to achieve the identification of the parameters $\\{\gamma_{l}\\}$.

The final effect function is calculated as follows:

$\hat{\beta}(t)=\boldsymbol{b}^{T}(t)\hat{\boldsymbol{\gamma}}=\sum_{l=1}^{L}b_{l}(t)\hat{\gamma}_{l}$
(10)

where the estimates $\hat{\boldsymbol{\gamma}}$ represent the MPCMR estimates
with transformed principal components. The pointwise standard error is denoted
as
$\hat{\sigma}(t)=\sqrt{\boldsymbol{b}^{T}(t)\hat{\Sigma}\boldsymbol{b}(t)}$.
The MVMR model with the general transformed principal components is
demonstrated in Figure 3.

Figure 3: The diagram of Mendelian randomization under the MPCMR model with
transformed principal components. $Z$ represents (high-dimensional) genetic
instruments, ${\xi_{1},\ldots,\xi_{K}}$ are the original $K$ principal
components for the individual exposure curve,
${\xi^{\ast}_{1},\ldots,\xi^{\ast}_{L}}$ are the transformed $L$ principal
components, $\vec{U}=\\{U(t);0\leq t\leq T\\}$ represents unmeasured
confounders, and $Y$ is the outcome.
$<\phi_{k},b_{l}>=\int\phi_{k}(t)b_{l}(t)dt$ for any $k$ and $l$, where
$\phi_{k}(t)$ is the $k$-th eigenfunction and $b_{l}(t)$ is the $l$-th known
basis function. ${\gamma_{1},\ldots,\gamma_{L}}$ are the parameters of
interest.

Note that one special example is to let
$\boldsymbol{b}(t)=\boldsymbol{\phi}(t)$; that is, the eigenfunction is also
the basis function. According to the eigenfunction property,
$\boldsymbol{B}=\boldsymbol{I}$, therefore the transformed principal
components are the principal components themselves. In addition to the
eigenfunctions, another classical choice of the basis functions are the
polynomial basis functions $b_{l}(t)=t^{l-1}$ for $l=1,2,\ldots$.

The approach we consider could be implemented in a two-sample setting, where
genetic variant associations with the exposure are measured in a non-
overlapping, but representative, sample to genetic variant associations with
the outcome, and with summary-level data, given that summary statistics
relating to the principal components are derived. The details for these
settings are provided in Supplementary Text S4-S5.

###  3.3 Time-varying IV strength analysis

The usual methods for assessing instrument relevance and validity can be used
for our time-varying MR models. For example, we can assess the evidence on
whether all the instruments used are valid using heterogeneity or over-
identification tests [12]. Similarly, we can evaluate instrument strength to
determine if the instrument can support meaningful and reliable time-varying
results. The IV strength or validity assessment for time-varying MR is similar
to those for MVMR, where the principal components are treated as the
exposures. Details are given in Supplementary Text S6-S7.

Instrument strength for the pseudo-exposure $\xi$ is closely related to
instrument strength for the time-varying exposure $X(t)$. Instead of measuring
the strength of genetic associations with the exposure at multiple time
points, the time-varying association can be evaluated in a more efficient way
using the FPCA approximation:

$\alpha(t)=\frac{cov(G,X(t))}{var(G)}\approx\frac{cov(G,\sum_{k=1}^{K}\xi_{k}\phi_{k}(t))}{var(G)}=\sum_{k=1}^{K}\frac{cov(G,\xi_{k})}{var(G)}\phi_{k}(t)=\sum_{k=1}^{K}{\alpha_{k}}\phi_{k}(t)$
(11)

where $\alpha_{k}$ represents the genetic association with the $k$-th
principal component. Therefore, genetic associations with the time-varying
exposure can be expressed as a weighted average of genetic associations with
the principal components. Instrument strength in multivariable MR is measured
by conditional F-statistics [38] (or Cochran’s Q statistic [36], both embedded
in our package TVMR).

###  3.4 Identification-robust MR analysis

If the time-varying genetic associations with the exposure are be similar
(i.e., approximately proportional) to each other, the genetic associations
with the principal components would be approximately proportional. This would
lead to the genetic variants being conditionally weak instruments. If we
increase the number of principal components used to model the exposure,
instrument strength is likely to get weaker. Therefore, we consider
identification-robust inference, to ensure that the probability of falsely
rejecting the null hypothesis can be well-controlled regardless of whether the
instrument is weak or the model is ill-identified [3]. Identification-robust
methods do not assume the exposure effects are point-identifiable and hence do
not return point estimates. Importantly, they return valid confidence
intervals even under weak instrument settings. In contrast, inferential
methods based on IVW or GMM point estimates can be biased under weakly
identified exposure effects. There are many statistics proposed for
identification-robust inference (see [3]). As a popular statistic that can
also work for over-identified models, Kleibergen’s Lagrange multiplier (LM)
statistic [17] is used in this paper. We provide a derivation to calculate and
perform inferences on this statistic based on summarized genetic association
data.

Following usual GMM arguments [17], we have the joint asymptotic distribution

$\sqrt{n}\left(\begin{array}[]{c}\hat{\boldsymbol{g}}(\boldsymbol{\beta}^{\ast})\\\
\textrm{vec}(\hat{\boldsymbol{G}})\end{array}\right)\overset{D}{\to}\mathcal{N}\left(\left(\begin{array}[]{c}\boldsymbol{0}\\\
\boldsymbol{\mu}\end{array}\right),\left(\begin{array}[]{cc}\Omega(\boldsymbol{\beta}^{\ast})&\Delta\\\
\Delta^{T}&\Sigma_{G}\end{array}\right)\right)$ (12)

where $\text{vec}(\hat{\boldsymbol{G}})$ is the $JL$-vector obtained by
stacking the columns of $\hat{\boldsymbol{G}}$. We let
$\Delta\in\mathbb{R}^{J\times JL}$ denote the covariance matrix of
$\hat{\boldsymbol{g}}(\boldsymbol{\beta}^{\ast})$, and let
$\Sigma_{G}\in\mathbb{R}^{JL\times JL}$ denote the the covariance matrix of
$\textrm{vec}(\hat{\boldsymbol{G}})$. The vector
$\sqrt{n}\textrm{vec}(\hat{\boldsymbol{G}})$ has a fixed mean vector
$\boldsymbol{\mu}$ due to the weak instrument assumption that the covariance
between any exposure (principal component) and instruments is decreasing at
the same rate as the sampling error $1/\sqrt{n}$. Under the null hypothesis
$H_{0}:\boldsymbol{\beta}^{\ast}=\boldsymbol{\beta}_{0}$, the weak-instrument-
robust approach aims to construct a test statistic (LM statistic) based on a
sufficient statistic for $\hat{\boldsymbol{G}}$ that is asymptotically
uncorrelated with $\hat{\boldsymbol{g}}(\boldsymbol{\beta}_{0})$ under
$H_{0}$. The calculation of LM statistic under $H_{0}$ involves three steps
(for more details, see the supplementary text S8):

  * 1.

We construct a consistent estimator of $\Delta$. The $k$-th column of
$\hat{\boldsymbol{G}}$ is
$\hat{\boldsymbol{G}}_{k}=\hat{\mathbb{E}}_{n}(-\boldsymbol{Z}\xi_{k})$.
Denote $\Delta=(\Delta_{1}\,\Delta_{2}\,\ldots\,\Delta_{K})$ where
$\Delta_{k}:=cov(\sqrt{n}\hat{\boldsymbol{g}}(\boldsymbol{\beta}_{0}),\sqrt{n}\hat{\boldsymbol{G}}_{k})$.
We can construct consistent estimators
$\hat{\Delta}_{k}=\hat{\mathbb{E}}_{n}(-\boldsymbol{Z}\xi_{k}(Y-\boldsymbol{\xi}^{T}\boldsymbol{\beta}_{0})\boldsymbol{Z}^{T})$
and $\hat{\Delta}=(\hat{\Delta}_{1}\,\ldots\,\hat{\Delta}_{K})$.

  * 2.

Then, for each principal component k, let
$\hat{D}_{k}(\boldsymbol{\beta}_{0})=\hat{\boldsymbol{G}}_{k}-\hat{\Delta}_{k}(\boldsymbol{\beta}_{0})^{T}\hat{\Omega}(\boldsymbol{\beta}_{0})^{-1}\hat{\boldsymbol{g}}(\boldsymbol{\beta}_{0})$.
A sufficient statistic of $\hat{\boldsymbol{G}}$ is
$\hat{\boldsymbol{D}}(\boldsymbol{\beta}_{0})=(\hat{D}_{1}(\boldsymbol{\beta}_{0})\,\ldots\,\hat{D}_{K}(\boldsymbol{\beta}_{0}))$;
it can be shown that each principal component $k$,
$\hat{D}_{k}(\boldsymbol{\beta}_{0})$ is asymptotically uncorrelated with
$\hat{\boldsymbol{g}}(\boldsymbol{\beta}_{0})$.

  * 3.

Finally, the LM statistic is

$LM(\boldsymbol{\beta}_{0})=[\hat{\Omega}(\boldsymbol{\beta})_{0}^{-\frac{1}{2}}\sqrt{n}\hat{\boldsymbol{g}}(\boldsymbol{\beta}_{0})]^{T}\hat{\boldsymbol{P}}(\boldsymbol{\beta}_{0})[\hat{\Omega}(\boldsymbol{\beta}_{0})^{-\frac{1}{2}}\sqrt{n}\hat{\boldsymbol{g}}(\boldsymbol{\beta}_{0})]$
(13)

where
$\hat{\boldsymbol{P}}(\boldsymbol{\beta}_{0})=[\hat{\Omega}(\boldsymbol{\beta})_{0}^{-\frac{1}{2}}\hat{\boldsymbol{D}}(\boldsymbol{\beta}_{0})][\hat{\boldsymbol{D}}(\boldsymbol{\beta}_{0})^{T}\hat{\Omega}(\boldsymbol{\beta}_{0})^{-1}\hat{\boldsymbol{D}}(\boldsymbol{\beta}_{0})]^{-1}[\hat{\Omega}(\boldsymbol{\beta})_{0}^{-\frac{1}{2}}\hat{\boldsymbol{D}}(\boldsymbol{\beta}_{0})]^{T}$

Under $H_{0}:\boldsymbol{\beta}^{\ast}=\boldsymbol{\beta}_{0}$ and some
regularity conditions, Theorem 1 of [17] shows that
$LM(\boldsymbol{\beta}_{0})\overset{D}{\to}\chi^{2}_{K}$. Thus, confidence
regions for $\beta^{\star}$ can be obtained by inverting the LM test. In
particular, one can follow these steps: (i) Assign $m$ increasing points for
each parameter, resulting in $m^{K}$ candidate values for
$\boldsymbol{\beta}_{0}$; (ii) Calculate the LM statistic over each candidate
value; (iii) Collect the parameter values for $\boldsymbol{\beta}^{\star}$
that are not rejected by an LM test. Candidates within the confidence region
for $\boldsymbol{\beta}^{\star}$ can be further used to obtain the confidence
interval for the effect function $\beta(t)$.

##  4 Simulation

###  4.1 Simulation design

The aim of our simulation study is to investigate the ability of the MPCMR
approach to estimate different forms of time-varying exposure effects. We
consider different scenarios for the structural exposure and outcome models.
Let the time region of interest be $[0,50]$. For the exposure, we consider the
following model:

$X(t)=\underbrace{\sum_{j=1}^{30}\alpha_{j}(t)G_{j}}_{\textrm{gene
score}}+\underbrace{[U_{0}+U(t)]}_{\textrm{confounding
effect}}+\epsilon_{X}(t)\qquad t\in[0,50]$ (14)

where $G_{j}$ represents the $j$-th genetic variant and $\alpha_{j}(t)$ is the
time-varying genetic effect on the exposure. Genetic variants were generated
as $G_{j}\overset{\textrm{i.i.d}}{\sim}B(2,0.3)$. The terms $U(t)$ and
$\epsilon_{X}(t)$ were set to be independent Wiener processes with variance
equal to $1$ at $t=50$, and $U_{0}\sim\mathcal{N}(0,1^{2})$ is a time-fixed
confounder. We consider three functions for the instrument–exposure
relationships, a sinusoidal function (A) and two linear functions (B and C):

  1. A.

$\alpha_{j}(t)=0.05\sin{(a_{j}t)}+b_{j}$ where
$a_{j}\overset{\textrm{i.i.d}}{\sim}Unif(-0.1,0.1)$ and
$b_{j}\overset{\textrm{i.i.d}}{\sim}Unif(-0.1,0.1)$

  2. B.

$\alpha_{j}(t)=a_{j}+b_{j}t$ where
$a_{j}\overset{\textrm{i.i.d}}{\sim}Unif(-0.1,0.1)$ and
$b_{j}\overset{\textrm{i.i.d}}{\sim}Unif(-0.004,0.004)$

  3. C.

$\alpha_{j}(t)=a_{j}+b_{j}t$ where
$a_{j}\overset{\textrm{i.i.d}}{\sim}Unif(-0.1,0.1)$ and
$b_{j}\overset{\textrm{i.i.d}}{\sim}Unif(-0.01,0.01)$

Figure 4: Examples of the individual variables for the exposure, gene score,
and confounding curves over time under three different instrument-exposure
models A-C.

The instrument strength for the exposure at any single timepoint is controlled
for scenarios A-B such that the average proportion of variation explained by
genetic variants is approximately $0.05$. We also control the genetic time-
varying effects so that different scenarios have different average instrument
strength measurements in the time-varying MR context. This will help to
investigate the performance of our method under different IV strengths. Figure
4 shows examples of the time-varying variables under different scenarios. We
selected $K=2$ principal components that explained at least 95% of exposure
variation. The corresponding first two eigenfunctions are of an approximately
linear form. We consider two sets of basis functions for estimating the time-
varying effect of the exposure on the outcome: the eigenfunctions, and
polynomial functions of degree 1 (i.e. linear functions). Fitting more complex
exposure–outcome models would require more principal components to be
estimated. However, this would lead to weaker instruments. We return to
discuss the trade-off between model complexity and instrument strength in
Section 6.

We consider the outcome model

$Y=\int_{0}^{50}\beta(t)X(t)dt+10[U_{0}+U(50)]+\epsilon_{Y}(50)$ (15)

We consider the following six exposure-outcome models with different causal
shapes $\beta(t)$, denoted by Scenario 1-6:

  1. 1.

Lifetime null effect: $\beta(t)\equiv 0$

  2. 2.

Lifetime constant effect: $\beta(t)\equiv 0.1$

  3. 3.

Lifetime increasing effect: $\beta(t)=0.02t$

  4. 4.

Lifetime sign-change effect: $\beta(t)=0.5-0.02t$

  5. 5.

Early-age-only effect: $\beta(t)=0.05(-t+20)I\\{t<20\\}$

  6. 6.

Later-age-only effect: $\beta(t)=0.05(t-30)I\\{t>30\\}$

Scenario 1 is the joint sharp causal null case [42]. Scenario 2 is the
constant effect case, which means the intervention of the exposure at any
timepoint will have the same effect on the outcome. Scenario 3 indicates that
the exposure has an increasingly stronger effect on the outcome over age,
which means the later the intervention of exposure is given, the greater the
effect on the outcome. Scenario 4 indicates that early-age exposure has a
positive effect but later-age exposure has a negative effect. Scenarios 5 and
6 are threshold effect cases where there exists a certain period under which
intervening on the exposure influences the outcome. Note that the outcome
models 1-4 are of a linear form, which can be well captured by the polynomial
basis function as well as the eigenfunctions (since the first two
eigenfunctions are approximately linear). However, the non-smooth effect
functions in scenarios 5-6 cannot be well approximated by the polynomials or
eigenfunctions. Therefore, our approach is unlikely to perform well under
Scenarios 5 and 6. Instead, one should consider other appropriate basis
functions according to domain knowledge regarding the underlying shape of the
effect function, or consider data-driven approaches to select relevant basis
functions.

We consider $n=10,000$ individuals in the simulation and randomly select 10
time points for each individual from the region $[0,50]$ to mimic a sparse
measurement scenario, which is common in practice. We first calculate
functional principal components using the widely-used package ‘fdapace’ [48].
We then fit the time-varying model using three different strategies:

  1. $\bullet$

Association: naive association fitting for the principal components and the
outcome without using instruments;

  2. $\bullet$

MPCMR(eigenfunction): MPCMR where the basis functions are the eigenfunctions;

  3. $\bullet$

MPCMR(polynomial): MPCMR where the basis functions are linear functions.

We use the weak instrument-robust LM statistic to calculate confidence
intervals for the fitted effect function for strategies (ii) and (iii). When
building the confidence intervals with the LM statistic, for each parameter
$\beta^{\ast}_{k}$, we search over the parameter regions
$[\hat{\beta}^{\ast}_{k}-4s.e.(\hat{\beta}^{\ast}_{k}),\hat{\beta}^{\ast}_{k}+4s.e.(\hat{\beta}^{\ast}_{k})]$,
which is sufficiently large to cover non-rejected regions in our simulation
setting.

We simulate and fit the data $1,000$ times independently for each scenario and
calculate the Mean Squared Error (MSE) of the GMM point estimates and coverage
for the effect function rate at $4$ time points ($t=10,20,30,40$). As is
common in functional data analysis, we avoid comparisons over time points near
the boundaries (e.g., $t=0,50$), where the curve estimates tend to have much
greater errors [28, Page 68]. To measure instrument strength, we calculate the
conditional F-statistic relating to each principal component.

###  4.2 Simulation results

We first discuss the exposure-outcome model scenario 3 (the lifetime
increasing effect) as a worked example. The time-varying MR results for
different instrument-exposure models and estimation strategies for one single
simulated dataset are shown in Figure 5. For the naive association approach
that does not use instruments, a persistent bias is observed over time due to
confounding, as the principal components are correlated with the confounders.
With MPCMR, the fitted effect function curve, whether using eigenfunctions or
linear polynomials as basis functions, closely resembles the true effect
function. For MPCMR fitting, compared with the conventional confidence
intervals based on GMM point estimates (denoted by the dashed black curves),
LM-based confidence intervals (denoted by the dashed grey curves) are more
conservative. The true effect function is covered by LM-based confidence
intervals at all timepoints in this simulated example. MPCMR exhibits good
performance in this example, irrespective of the instrument-exposure model
(sinusoidal or linear). Note that MPCMR using the two different basis
functions produces similar effect function estimates. This is because the
eigenfunctions used (the first two eigenfunctions) are both approximately
linear (see Supplementary Figure A1), resulting in similar linear results for
the two MPCMR strategies.

Now we examine the performance of different estimating strategies in time-
varying MR for all the exposure and outcome model scenarios considered.
Results are provided in Figure 6 and Table 1. In all the considered
timepoints, MPCMR (using eigenfunction and polynomial) shows better
performance in terms of MSE, and coverage rates compared to naive estimation
without using instruments (i.e. the association fitting) for outcome model
scenarios 1-4, where the true effect function is linear and can be well
approximated by eigenfunctions or polynomials. The coverage rate using LM
statistics is at an appropriate level (approximately $0.95$) in most cases.
The three exposure models (A-C) exhibit increasing instrument strength,
reflected by the increasing conditional F values for each principal component.
With an increase in IV strength, MPCMR usually performs better in terms of
smaller MSE. When the effect function is linear, MPCMR estimation using the
correct basis function (e.g., the polynomial basis) has better MSE performance
than MPCMR using eigenfunctions as the basis function, which are determined by
the exposure data rather than user-defined.

When the true effect function cannot be well approximated by the linear
combination of basis functions (corresponding to outcome models 5-6), our
effect function estimates have larger biases at some timepoints (also see
Supplementary Figure A2). Therefore, in addition to the IV assumptions, one
should also consider the functional form of the underlying effect function and
choose appropriate basis functions. Adhering to the aphorism that all models
are wrong but some are useful [6], the linear basis function may still offer
valuable insights into determining the overall trend of time-varying effects
and, consequently, the relative importance between the exposure over early and
late ages on the outcome. Taking outcome scenario 5 and the left panel of
Supplementary Figure A2 as an example, the negative slope of the fitted
function supports the notion that the exposure during early ages is more
likely to have a more positive effect than exposure during later ages.
However, as the models are misspecified, there is considerable bias in
estimation in the threshold outcome scenarios.

Figure 5: Time-varying MR fitting results using different estimation
strategies under different instrument-exposure models under one single
simulated example. The rows correspond to the instrument-exposure models A, B,
and C, respectively. The columns correspond to the estimation strategies:
association fitting without using instruments, MPCMR fitting where the basis
function is the eigenfunctions (therefore the exposure is not transformed),
and MPCMR fitting where the basis function is the polynomial functions of
degree 1 (therefore with a transformed exposure). The black solid and dashed
curves represent the estimated effect function results and their 95%
confidence intervals from the GMM method. The grey dashed curves represent the
confidence intervals built by the Lagrange multiplier statistic. The
underlying time-varying effect is the lifetime increasing effect (Scenario 3),
which is represented by a blue line.  Figure 6: Mean time-varying predicted
effect values using different estimation strategies under different
instrument-exposure models with 1,000 simulations at 6 timepoints
($t=0,10,20,30,40,50$). The columns correspond to the instrument-exposure
models A, B, and C, respectively. The rows correspond to the exposure-outcome
models 1-6. Association: the association fitting without using IVs.
MPCMR(eigenfunction): MPCMR fitting where the basis functions are the
eigenfunctions. MPCMR(polynomial): MPCMR fitting where the basis functions are
linear polynomials.

| Scenario A | Scenario B | Scenario C
---|---|---|---
| 8.043 (PC1) 7.284 (PC2) | 16.414 (PC1) 13.861 (PC2) | 25.918 (PC1) 20.900 (PC2)
| A1 | A2 | A3 | A4 | A5 | A6 | B1 | B2 | B3 | B4 | B5 | B6 | C1 | C2 | C3 | C4 | C5 | C6
$t=10$ | MSE ($10^{-2}$) | Association | 1.149 | 1.166 | 1.047 | 1.361 | 0.176 | 0.022 | 1.694 | 1.722 | 1.564 | 1.990 | 0.479 | 0.140 | 4.261 | 4.325 | 4.212 | 4.648 | 2.124 | 1.313
MPCMR(eigenfunction) | 0.456 | 0.468 | 0.757 | 0.374 | 3.002 | 4.630 | 0.321 | 0.345 | 0.259 | 0.928 | 0.464 | 0.857 | 0.327 | 0.364 | 0.206 | 0.931 | 0.460 | 0.551
MPCMR(polynomial) | 0.445 | 0.442 | 0.548 | 0.378 | 3.006 | 4.140 | 0.310 | 0.345 | 0.261 | 0.763 | 0.405 | 0.696 | 0.324 | 0.375 | 0.257 | 0.832 | 0.433 | 0.459
Coverage (%) | Association | 0.000 | 0.000 | 0.000 | 0.000 | 3.300 | 55.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.600 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000
MPCMR(eigenfunction) | 98.400 | 98.400 | 97.300 | 96.800 | 70.700 | 51.400 | 98.700 | 98.500 | 96.000 | 91.500 | 97.000 | 66.800 | 98.300 | 98.200 | 95.800 | 91.000 | 97.200 | 56.600
MPCMR(polynomial) | 95.300 | 95.500 | 96.100 | 94.500 | 63.200 | 48.100 | 97.100 | 97.000 | 97.100 | 89.000 | 95.400 | 69.900 | 98.000 | 97.800 | 97.600 | 90.600 | 97.000 | 62.200
$t=20$ | MSE ($10^{-2}$) | Association | 2.192 | 2.433 | 1.764 | 4.132 | 22.229 | 4.382 | 2.543 | 2.794 | 2.166 | 4.430 | 23.160 | 5.104 | 3.765 | 4.084 | 3.485 | 5.845 | 26.465 | 6.953
MPCMR(eigenfunction) | 0.117 | 0.110 | 0.256 | 0.403 | 8.311 | 0.148 | 0.143 | 0.194 | 0.110 | 1.058 | 15.517 | 0.820 | 0.150 | 0.214 | 0.091 | 1.028 | 15.586 | 1.000
MPCMR(polynomial) | 0.067 | 0.066 | 0.093 | 0.066 | 6.147 | 0.404 | 0.092 | 0.101 | 0.070 | 0.242 | 10.947 | 1.293 | 0.104 | 0.117 | 0.072 | 0.274 | 11.277 | 1.430
Coverage (%) | Association | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000
MPCMR(eigenfunction) | 98.100 | 97.600 | 95.000 | 76.900 | 0.200 | 96.200 | 98.800 | 98.700 | 94.600 | 60.600 | 0.000 | 74.500 | 98.500 | 98.100 | 95.500 | 59.300 | 0.000 | 62.200
MPCMR(polynomial) | 95.200 | 95.400 | 95.600 | 92.900 | 0.000 | 55.100 | 96.700 | 96.900 | 96.100 | 86.700 | 0.000 | 17.400 | 97.300 | 97.300 | 96.500 | 89.100 | 0.000 | 15.100
$t=30$ | MSE ($10^{-2}$) | Association | 3.576 | 3.757 | 4.032 | 4.064 | 10.025 | 22.024 | 3.325 | 3.487 | 3.795 | 3.702 | 9.437 | 21.553 | 2.432 | 2.575 | 2.727 | 2.888 | 8.046 | 18.908
MPCMR(eigenfunction) | 0.104 | 0.140 | 0.229 | 0.218 | 3.096 | 10.370 | 0.019 | 0.032 | 0.103 | 0.057 | 1.901 | 9.103 | 0.011 | 0.022 | 0.058 | 0.061 | 1.972 | 8.668
MPCMR(polynomial) | 0.081 | 0.083 | 0.053 | 0.145 | 2.465 | 9.483 | 0.014 | 0.014 | 0.011 | 0.024 | 1.501 | 8.291 | 0.007 | 0.007 | 0.006 | 0.020 | 1.506 | 8.175
Coverage (%) | Association | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000
MPCMR(eigenfunction) | 97.800 | 97.100 | 94.000 | 88.100 | 1.000 | 0.000 | 98.900 | 97.300 | 71.800 | 83.200 | 0.000 | 0.000 | 97.900 | 95.600 | 66.600 | 70.000 | 0.000 | 0.000
MPCMR(polynomial) | 95.000 | 94.700 | 93.600 | 93.100 | 3.000 | 0.000 | 57.100 | 56.300 | 47.600 | 56.600 | 0.000 | 0.000 | 85.800 | 85.900 | 68.600 | 65.100 | 0.000 | 0.000
$t=40$ | MSE ($10^{-2}$) | Association | 4.928 | 4.808 | 5.370 | 3.962 | 2.300 | 4.316 | 3.924 | 3.818 | 4.277 | 3.107 | 1.597 | 3.286 | 1.275 | 1.213 | 1.391 | 0.898 | 0.199 | 0.887
MPCMR(eigenfunction) | 0.687 | 0.696 | 1.173 | 0.442 | 0.728 | 1.710 | 0.098 | 0.100 | 0.229 | 0.418 | 1.495 | 0.106 | 0.048 | 0.057 | 0.060 | 0.367 | 1.505 | 0.060
MPCMR(polynomial) | 0.489 | 0.491 | 0.426 | 0.613 | 0.730 | 0.645 | 0.077 | 0.082 | 0.085 | 0.111 | 0.813 | 0.193 | 0.035 | 0.044 | 0.058 | 0.070 | 0.833 | 0.217
Coverage (%) | Association | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 2.800 | 0.000
MPCMR(eigenfunction) | 98.900 | 98.900 | 97.000 | 97.000 | 96.800 | 93.400 | 98.800 | 98.700 | 92.500 | 77.500 | 19.700 | 98.900 | 98.900 | 98.300 | 83.000 | 58.000 | 1.600 | 97.500
MPCMR(polynomial) | 96.400 | 96.400 | 96.500 | 96.500 | 93.200 | 94.700 | 97.300 | 97.400 | 97.500 | 95.100 | 32.600 | 90.400 | 97.900 | 98.000 | 96.800 | 95.000 | 3.300 | 65.000

Table 1: Time-varying MR fitting results using three estimation strategies
with different combinations of the instrument-exposure models (A-C) and the
exposure-outcome models (1-6). The results include the mean squared error
(MSE) based on the GMM point estimates and the average coverage rate at the
four timepoints ($t=10,20,30,40$) with 1,000 independent simulations. For each
instrument-exposure model, the average conditional F statistics for each
principal component were recorded. PC1 and PC2 represent the first and second
principal components, respectively. Association represents association fitting
without using IVs. MPCMR(eigenfunction) represents MPCMR fitting where the
basis functions are the eigenfunctions. MPCMR(polynomial) represents MPCMR
fitting where the basis functions are linear polynomials. The coverage
calculations for MPCMR are based on the identification-robust IV statistics.

##  5 Illustrative example: Systolic Blood Pressure on Urea Level

We give an example of applying our time-varying MR method to a real data
application to illustrate how our method works. Urea is a waste product
removed from the blood by the kidneys. Chronic high blood pressure can impair
kidney function, which may affect urea levels.

We analyzed data from UK Biobank, a population-based cohort study of UK
residents ages 40-69 at recruitment [41], who passed quality control checks as
previously described [5]. We investigate the effect of systolic blood pressure
(SBP) on urea level using 258 genetic variants as instruments. These genetic
variants were previously shown to be associated with a blood pressure variable
in a genome-wide association study that did not include UK Biobank
participants [11], thus avoiding bias due to winner’s curse [44].

We first investigate the effect of SBP on urea level in a standard univariable
MR model ignoring time, where both SBP and urea level are treated as time-
invariant. The SBP value measured in mmHg at recruitment is used as the
exposure, and urea measured in mmol/L at recruitment is used as the outcome.
Univariable analyses are conducted in 349,883 participants, with blood
pressure and urea measured at study baseline. Supplementary Figure A3 shows
the scatter plot of the genetic associations with the exposure and the outcome
using the 258 SNPs. The univariable MR effect estimate is $0.008$ (95% CI:
$0.004,0.012$) mmol/L per $1$ mmHg higher genetically-predicted SBP. This
suggests that higher genetically-predicted SBP levels are associated with
higher urea levels. The estimate may be interpreted as a cumulative effect,
under the assumption that the genetic associations with the exposure are
approximately constant [18]. However, a positive cumulative effect cannot
provide insight on any time-varying effects. Therefore, we use the MPCMR
method to investigate the possible effect pattern over time.

For the time-varying analysis, we use SBP measurements from primary care
records on 158,310 individuals, and their urea measurements taken at baseline.
To ensure there was sufficient data to obtain precise estimates at all
timepoints, we considered blood pressure measurements from age 50 to 75 years.
Participants without a SBP measurement in this time window were dropped from
analyses. To reduce computational complexity, we rounded all measurements to
the nearest tenth of a year.

The data used for MPCMR fitting contains 101,279 individuals. We calculated
the eigenfunctions and individual principal components from FPCA, using the
package _fdapace_. The eigenfunction information is given in Supplementary
Figure A4. The first four principal components have a cumulative fraction of
variance explained larger than $95\%$, with the first, second, third, and
fourth components explaining $65.9\%$, $18.5\%$, $8.2\%$, and $4.7\%$ of the
variance in exposure measurements. We use the first two eigenfunctions and the
linear polynomials as the basis functions in MPCMR fitting to mitigate bias
from weak instruments.

The functional fitting results are given in Figure 7. The estimated effect
(represented by the solid black curve) is decreasing over time and significant
and positive over some time points. Both eigenfunction and linear polynomials
fitting with the GMM 95% confidence interval (the black dashed curves) support
that the early-age effect of SBP is positive and significant, while the later-
age effect is not significant. This indicates that SBP at the beginning of the
time window is more likely to have a more positive effect than SBP towards the
end of the time window. However, the results still indicate much uncertainty,
and the estimated effect is significant for only a short period around 60
years as indicated by the identification-robust 95% confidence intervals (the
grey dashed curves). The high uncertainty is due to the weak-instrument
problem as there is a strong pattern of linear dependence in the genetic
associations with the first and second principal components. The conditional F
statistics are $9.46$ and $1.52$ for the two principal components. Overall, we
cannot made any strong conclusions about the time-varying nature of the effect
of SBP on urea levels.

Figure 7: Application results using the first two eigenfunctions (left panel)
and the linear polynomials (right panel) as the basis functions. The black
curve represents the fitted time-varying effect from the MPCMR method, while
the dashed lines represent 95% confidence intervals derived from GMM (black)
and Lagrange multipliers (grey). The effect unit is mmol/L per $1$ mmHg higher
genetically-predicted SBP.

##  6 Discussion

In this paper, we have introduced a novel method for Mendelian randomization
to study the time-varying effects of a modifiable exposure on the outcome. Our
method is based on continuous-time modelling with a functional effect
objective, so it should not suffer from violation of the IV exclusion
restriction assumption caused by the model misspecification arising from using
exposure measurements at limited measured timepoints. A continuous-time model
has previously been proposed for time-varying MR studies [9, 18] but there was
no feasible estimation method proposed. We have proposed a general framework
for dealing with such problem relying on functional dimensional reduction
ideas. Our method initially constructs a multivariable MR model using
functional principal component analysis summarizing time-varying exposure
information. We then back out time-varying effect estimation and inference
based on the MVMR model estimates. We considered weak-instrument-robust
methods so that we can obtain reliable results regardless of instrument
strength. We demonstrated with a variety of simulation scenarios that our
method can accurately estimate time-varying effect functions, given that the
correct form of the basis function is used.

To the best of our knowledge, all the current methods proposed in time-varying
MR studies are time-discrete. For example, a previous paper considered effect
estimation using multiple exposure values, but it only estimated a constant
causal effect [9]. Other past papers assumed the effect function can be
correctly expressed as a constant effect at limited timepoints [39, 37]. All
of those additional assumptions can lead to misspecification problems and
yield misleading MR results [18, 45]. By contrast, we treat the time-varying
effect function as the objective of interest since time is continuous, and so
our approach can provides a more general way for conducting time-varying
analyses.

In order to estimate the effect function, we perform dimension reduction to
transform the exposure trajectories (or equivalently the genetic associations
with the exposure) into nonfunctional variables. We employed FPCA, a data-
driven method that borrows information across individuals and avoids making
strong parametric model assumptions for the exposure trajectory or instrument-
exposure associations. Our approach allows for sparse and different exposure
measurement time points for each individual, which is common in practice. We
adopted identification-robust inference and have demonstrated that valid
conclusions can be obtained by our method for linear time-varying functions.
It is important to consider weak instrument-robust inference in time-varying
MR analyses, since it may be difficult to find genetic predictors of the
exposure at multiple time points with appropriate linear independence.

Although our MPCMR approach allows for a flexible parametric form to learn the
effect function, we have presented simulation scenarios using eigenfunction
and linear basis functions, which are adequate for estimating linear effect
functions, but are unable to estimate more complex functions. This choice is
based on practical considerations as FPCA may only return a relatively small
number of principal components explaining the exposure trajectory. Increasing
the number of principal components used for MPCMR can worsen the weak-
instrument problem, as the added principal components explain a lower
proportion of exposure variation and so could be less associated with the
instruments. This presents a trade-off in MPCMR estimation: if one desires a
more flexible estimate of the time-varying function, they may suffer more from
weak instrument issues. While if one uses fewer PCs to enable model fitting
with strong IVs, this requires use of a simpler form for fitting the effect
function. In particular, if we only have two principal components, then we can
only estimate two parameter exposure–outcome models. This means that
estimating complex time-varying functions is likely to be infeasible in
practice. With the increasing availability of detailed time-stamped diagnoses
and longitudinal data from electronic health records in large biobanks, such
as UK Biobank [8] and FinnGen (https://www.finngen.fi/en), we believe time-
varying MR analysis has the potential to become more widely implemented in
future studies. However, this will require the existence of genetic variants
having effects on exposure trajectories that vary in time. This is a natural
limitation to the feasibility of any approach for time-varying MR studies.

There are some methodological limitations in this study. First, to extract
reliable principal component information, FPCA requires individual exposure
values measured at multiple time points, with each time point containing a
sufficient number of samples. This means the data should be longitudinal over
the time region of interest. Second, FPCA may not explain a sufficient
proportion of the variation in the exposure with the desired number of
principal components, and the remaining ignored principal components could
induce structural pleiotropy, which may be a concern for instrument validity.
In our example, the first two principal components explained around 85% of the
variability in the exposure trajectories, which is good from the perspective
of avoiding structural pleiotropy, but only allows the estimation of simple
exposure–outcome functions. Third, as in any application of MR, some
instruments may have direct effects on the outcome and may invalidate our
estimation, which requires the IV assumptions to hold. Additionally, using a
truncated, clipped, or incomplete time region may violate the IV assumptions,
even if all IVs are valid in the original complete time region. Fourth,
uncertainty in estimating the eigenfunctions is ignored in our approach.
Uncertainty can be reduced by increasing the sample size or the number of
measured time points. Additionally, we assume that the same eigenfunctions
hold for all individuals. Finally, although our method aims to be generally
applicable in time-varying MR studies, it still requires some assumptions to
ensure the model is correct. Future methodological work should consider the
relaxation of these assumptions.

In summary, we have presented a novel method with continuous-time modelling
that combines data-driven approaches and weak IV robust methods to obtain
reliable estimators for simple time-varying MR effects. Our method offers a
general framework for studying the effects of time-varying modifiable
exposures, but is limited by the availability of IVs with time-varying effects
on the exposure. Since most exposures of interest are time-varying, our method
has the potential to be an important element for future MR analyses. However,
there are intrinsic limitations on the level of detail for time-varying
effects that can learned in practice.

## Acknowledgements

For the purpose of open access, the authors have applied a Creative Commons
Attribution (CC BY) licence to any Author Accepted Manuscript version arising
from this submission.

## Code availability

The R package to apply the methods proposed in this manuscript is available at
https://github.com/HDTian/TVMR, along with R code that reproduces the
empirical results in this manuscript.

## References

  * [1] Odd O Aalen, Kjetil Røysland, Jon Michael Gran, and Bruno Ledergerber.  Causality, mediation and time: a dynamic viewpoint.  Journal of the Royal Statistical Society: Series A (Statistics in Society), 175(4):831–861, 2012.
  * [2] Odd Olai Aalen, Kjetil Røysland, Jon Michael Gran, Roger Kouyos, and Tanja Lange.  Can we believe the dags? a comment on the relationship between causal dags and mechanisms.  Statistical methods in medical research, 25(5):2294–2314, 2016.
  * [3] Isaiah Andrews, James H Stock, and Liyang Sun.  Weak instruments in instrumental variables regression: Theory and practice.  Annual Review of Economics, 11:727–753, 2019.
  * [4] Bertille Antoine and Eric Renault.  Efficient gmm with nearly-weak instruments.  The Econometrics Journal, 12:S135–S171, 2009.
  * [5] William J Astle, Heather Elding, Tao Jiang, Dave Allen, Dace Ruklisa, Alice L Mann, Daniel Mead, Heleen Bouman, Fernando Riveros-Mckay, Myrto A Kostadima, et al.  The allelic landscape of human blood cell trait variation and links to common complex disease.  Cell, 167(5):1415–1429, 2016.
  * [6] George EP Box.  Science and statistics.  Journal of the American Statistical Association, 71(356):791–799, 1976.
  * [7] Stephen Burgess and Simon G Thompson.  Mendelian randomization: methods for causal inference using genetic variants.  CRC Press, 2021.
  * [8] Clare Bycroft, Colin Freeman, Desislava Petkova, Gavin Band, Lloyd T Elliott, Kevin Sharp, Allan Motyer, Damjan Vukcevic, Olivier Delaneau, Jared O’Connell, et al.  The uk biobank resource with deep phenotyping and genomic data.  Nature, 562(7726):203–209, 2018.
  * [9] Ying Cao, Suja S Rajan, and Peng Wei.  Mendelian randomization analysis of a time-varying exposure for binary disease outcomes using functional data analysis methods.  Genetic epidemiology, 40(8):744–755, 2016.
  * [10] George Davey Smith and Shah Ebrahim.  ‘mendelian randomization’: can genetic epidemiology contribute to understanding environmental determinants of disease?  International journal of epidemiology, 32(1):1–22, 2003.
  * [11] Evangelos Evangelou, Helen R Warren, David Mosen-Ansorena, Borbala Mifsud, Raha Pazoki, He Gao, Georgios Ntritsos, Niki Dimou, Claudia P Cabrera, Ibrahim Karaman, et al.  Genetic analysis of over 1 million people identifies 535 new loci associated with blood pressure traits.  Nature genetics, 50(10):1412–1425, 2018.
  * [12] Lars Peter Hansen.  Large sample properties of generalized method of moments estimators.  Econometrica: Journal of the econometric society, pages 1029–1054, 1982.
  * [13] Lars Peter Hansen, John Heaton, and Amir Yaron.  Finite-sample properties of some alternative gmm estimators.  Journal of Business & Economic Statistics, 14(3):262–280, 1996.
  * [14] Miguel A Hernán and James M Robins.  Instruments for causal inference: an epidemiologist’s dream?  Epidemiology, pages 360–372, 2006.
  * [15] Joel L Horowitz.  Ill-posed inverse problems in economics.  Annu. Rev. Econ., 6(1):21–51, 2014.
  * [16] Guido W Imbens and Donald B Rubin.  Causal inference in statistics, social, and biomedical sciences.  Cambridge University Press, 2015.
  * [17] Frank Kleibergen.  Testing parameters in gmm without assuming that they are identified.  Econometrica, 73(4):1103–1123, 2005.
  * [18] Jeremy A Labrecque and Sonja A Swanson.  Interpretation and potential biases of mendelian randomization estimates with time-varying exposures.  American journal of epidemiology, 188(1):231–238, 2019.
  * [19] Jeremy A Labrecque and Sonja A Swanson.  Commentary: Mendelian randomization with multiple exposures: the importance of thinking about time.  International journal of epidemiology, 49(4):1158–1162, 2020.
  * [20] Jeremy A Labrecque and Sonja A Swanson.  Age-varying genetic associations and implications for bias in mendelian randomization.  medRxiv, 2021.
  * [21] Debbie A Lawlor, Roger M Harbord, Jonathan AC Sterne, Nic Timpson, and George Davey Smith.  Mendelian randomization: using genes as instruments for making causal inferences in epidemiology.  Statistics in medicine, 27(8):1133–1163, 2008.
  * [22] Haben Michael, Yifan Cui, Scott A Lorch, and Eric J Tchetgen Tchetgen.  Instrumental variable estimation of marginal structural mean models for time-varying treatment.  Journal of the American Statistical Association, pages 1–12, 2023.
  * [23] Tim T Morris, Jon Heron, Eleanor CM Sanderson, George Davey Smith, Vanessa Didelez, and Kate Tilling.  Interpretation of mendelian randomization using a single measure of an exposure that varies over time.  International Journal of Epidemiology, 51(6):1899–1909, 2022.
  * [24] Whitney K Newey and Daniel McFadden.  Large sample estimation and hypothesis testing.  Handbook of econometrics, 4:2111–2245, 1994.
  * [25] Whitney K Newey and James L Powell.  Instrumental variable estimation of nonparametric models.  Econometrica, 71(5):1565–1578, 2003.
  * [26] Laura Pazzagli, Marie Linder, Mingliang Zhang, Emese Vago, Paul Stang, David Myers, Morten Andersen, and Shahram Bahmanyar.  Methods for time-varying exposure related problems in pharmacoepidemiology: an overview.  Pharmacoepidemiology and drug safety, 27(2):148–160, 2018.
  * [27] Grace M Power, Eleanor Sanderson, Panagiota Pagoni, Abigail Fraser, Tim Morris, Claire Prince, Timothy M Frayling, Jon Heron, Tom G Richardson, Rebecca Richmond, et al.  Methodological approaches, challenges, and opportunities in the application of mendelian randomisation to lifecourse epidemiology: A systematic literature review.  European Journal of Epidemiology, pages 1–20, 2023.
  * [28] James Ramsay and B.W. Silverman.  Functional Data Analysis, Second Edition.  John Wiley & Sons, Ltd, 2005.
  * [29] Tom G Richardson, Genevieve M Leyden, and George Davey Smith.  Time-varying and tissue-dependent effects of adiposity on leptin levels: a mendelian randomization study.  medRxiv, pages 2022–11, 2022.
  * [30] Tom G Richardson, Eleanor Sanderson, Benjamin Elsworth, Kate Tilling, and George Davey Smith.  Use of genetic variation to separate the effects of early and later life adiposity on disease risk: mendelian randomisation study.  bmj, 369, 2020.
  * [31] James M Robins.  The analysis of randomized and non-randomized aids treatment trials using a new approach to causal inference in longitudinal studies.  Health service research methodology: a focus on AIDS, pages 113–159, 1989.
  * [32] James M Robins.  Correcting for non-compliance in randomized trials using structural nested mean models.  Communications in Statistics-Theory and methods, 23(8):2379–2412, 1994.
  * [33] James M Robins.  Marginal structural models versus structural nested models as tools for causal inference.  In Statistical models in epidemiology, the environment, and clinical trials, pages 95–133. Springer, 2000.
  * [34] Donald B Rubin.  Estimating causal effects of treatments in randomized and nonrandomized studies.  Journal of educational Psychology, 66(5):688, 1974.
  * [35] Donald B Rubin.  Randomization analysis of experimental data: The fisher randomization test comment.  Journal of the American statistical association, 75(371):591–593, 1980.
  * [36] Eleanor Sanderson, George Davey Smith, Frank Windmeijer, and Jack Bowden.  An examination of multivariable mendelian randomization in the single-sample and two-sample summary data settings.  International journal of epidemiology, 48(3):713–727, 2019.
  * [37] Eleanor Sanderson, Tom G Richardson, Tim T Morris, Kate Tilling, and George Davey Smith.  Estimation of causal effects of a time-varying exposure at multiple time points through multivariable mendelian randomization.  PLoS Genetics, 18(7):e1010290, 2022.
  * [38] Eleanor Sanderson and Frank Windmeijer.  A weak instrument f-test in linear iv models with multiple endogenous variables.  Journal of econometrics, 190(2):212–221, 2016.
  * [39] Joy Shi, Sonja A Swanson, Peter Kraft, Bernard Rosner, Immaculata De Vivo, and Miguel A Hernán.  Mendelian randomization with repeated measures of a time-varying exposure: an application of structural mean models.  Epidemiology, 33(1):84–94, 2021.
  * [40] George Davey Smith.  Randomised by (your) god: robust inference from an observational study design.  Journal of Epidemiology & Community Health, 60(5):382–388, 2006.
  * [41] C Sudlow, J Gallacher, N Allen, V Beral, P Burton, J Danesh, et al.  UK Biobank: an open access resource for identifying the causes of a wide range of complex diseases of middle and old age.  PLOS Medicine, 12(3):e1001779, 2015.
  * [42] Sonja A Swanson, Jeremy Labrecque, and Miguel A Hernán.  Causal null hypotheses of sustained treatment strategies: What can be tested with an instrumental variable?  European journal of epidemiology, 33:723–728, 2018.
  * [43] Sonja A Swanson, Henning Tiemeier, M Arfan Ikram, and Miguel A Hernán.  Nature as a trialist? deconstructing the analogy between mendelian randomization and randomized trials.  Epidemiology (Cambridge, Mass.), 28(5):653, 2017.
  * [44] Amy E Taylor, Neil M Davies, Jennifer J Ware, Tyler VanderWeele, George Davey Smith, and Marcus R Munafò.  Mendelian randomization in health research: using appropriate genetic variants and avoiding biased estimates.  Economics & Human Biology, 13:99–106, 2014.
  * [45] Haodong Tian and Stephen Burgess.  Estimation of time-varying causal effects with multivariable mendelian randomization: some cautionary notes.  International journal of epidemiology, 52(3):846–857, 2023.
  * [46] Jane-Ling Wang, Jeng-Min Chiou, and Hans-Georg Müller.  Functional data analysis.  Annual Review of Statistics and its application, 3:257–295, 2016.
  * [47] Fang Yao, Hans-Georg Müller, and Jane-Ling Wang.  Functional data analysis for sparse longitudinal data.  Journal of the American statistical association, 100(470):577–590, 2005.
  * [48] Yidong Zhou, Satarupa Bhattacharjee, Cody Carroll, Yaqing Chen, Xiongtao Dai, Jianing Fan, Alvaro Gajardo, Pantelis Z. Hadjipantelis, Kyunghee Han, Hao Ji, Changbo Zhu, Hans-Georg Muller, and Jane-Ling Wang.  fdapace: Functional Data Analysis and Empirical Dynamics, 2022.  R package version 0.5.9.

## Supplementary Figures

Supplementary Figure A1: The eigenfunction results of functional principal
components in one simulation example under the three exposure model scenarios
A-C. The left, middle and right panel corresponds to the instrument-exposure
scenarios A, B and C, respectively. In all cases, the first and second
eigenfunctions, denoted by the red and green curves, are approximately linear
and cumulatively explain more than $95\%$ of the variations.

Supplementary Figure A2: The time-varying MR fitting results using
eigenfunctions as the basis functions under the exposure model scenario C.
Left: the outcome model scenario 5. Right: the outcome model scenario 6. The
true effect functions, which are expressed by the blue curves, cannot be
represented by the eigenfunctions in an additive linear way, thereby causing
bias in most timepoints. Supplementary Figure A3: Scatter plot of the
estimated genetic association with the exposure and outcome for 258 SNPs in
the univariable Mendelian randomization fitting. Supplementary Figure A4: The
eigenfunctions of FPCA for the SBP data over the age region $[50,75]$. PC
stands for the principal component, followed by the corresponding fraction
variance explained.

## Supplementary Texts

### Text S1: Basic identification equation

We illustrate the key idea of identifying the time-varying effect, which is a
functional objective, using instrumental variable (IV) methods. Assuming the
IV core assumptions are valid, and the structural outcome is

$Y_{T}=\int_{0}^{T}\beta_{T}(t)X(t)dt+g_{Y,T}(U,\epsilon_{Y})$ (S1)

which can be regarded as an infinitely dimensional multivariable Mendelian
randomization (MVMR) model, where the multiple exposures represent the time-
varying exposures at different time points. $U$ contains both the time-
invariant and time-varying confounders up to the time point $T$ (e.g.
$\\{U(t);0\leq t\leq T\\}$). Let’s assume there are $J$ independent genetic
variants denoted as ${G_{j};j=1,\ldots,J}$. The inverse variance weighted
(IVW) method involves two-stage fittings. In the second stage, we perform a
regression of the outcome on the instruments (with possible covariates for
adjustment). The estimated parameter for the $j$-th instrument follows an
asymptotic distribution (here, we ignore the uncertainty associated with the
standard errors):

$\hat{\theta}_{j}\overset{D}{\to}\mathcal{N}(\theta_{j},s.e.(\hat{\theta}_{j})^{2})$
(S2)

where
$\theta_{j}=cov(Y_{T},G_{j})/var(G_{j})=cov(\int_{0}^{T}\beta_{T}(t)X(t)dt,G_{j})/var(G_{j})=\int_{0}^{T}\beta_{T}(t)\alpha_{j}(t)dt$
and $\alpha_{j}(t):=cov(G_{j},X(t))/var(G_{j})$. It is easy to know
$cov(\hat{\theta}_{j_{1}},\hat{\theta}_{j_{2}})=0$ for any $j_{1}\neq j_{2}$.
Therefore, we have the following equation that enables the identification of
$\beta_{T}(t)$:
$\hat{\theta}_{j}\sim\mathcal{N}(\int_{0}^{T}\beta_{T}(t)\alpha_{j}(t)dt,s.e.(\hat{\theta}_{j})^{2})$,
which is the ideal fitting equation for IVW results. However, since
$\alpha_{j}(t)$ is unknown, it needs to be estimated from the data. This
estimation can be obtained through the first-stage IVW fitting, which involves
regressing $X(t)$ on the instruments (with possible covariates for
adjustment). If the uncertainty associated with the estimated functional
objective $\hat{\alpha}_{j}(t)$ is small, we can proceed with the following
fitting regression

$\
\hat{\theta}_{j}=\int_{0}^{T}\beta_{T}(t)\hat{\alpha}_{j}(t)dt+\epsilon_{j}\qquad\epsilon_{j}\sim\mathcal{N}(0,s.e.(\hat{\theta}_{j})^{2})\qquad
j=1,\ldots,J$ (S3)

which is the standard fitting equation for obtaining IVW results. Note that
this equation represents an ill-posed inverse problem, which means that
$\beta_{T}(t)$ may not be consistently estimated even if we have a consistent
estimator of $\theta_{j}$ and $\alpha_{j}(t)$ for $j=1,\ldots,J$. Therefore,
we cannot directly express the final estimated term (the effect function,
$\hat{\beta}_{T}(t)$, in time-varying MR) in a manner similar to univariable
MR or MVMR without making additional assumptions about $\beta_{T}(t)$.

### Text S2: Functional principal component analysis

We briefly introduce the functional principal component analysis (FPCA) and
its basic properties. For more details, see the review paper
wang2016functional. Assume the individual function curve
$X_{i}(t),i=1,\ldots,n$ is a smooth and square-integrable function, where
$X_{i}(t)$ is a functional objective like the exposure trajectory. The FPCA
treats $X_{i}(t)$ as a realization of the random function $X(t)$ over the
region of interest $[0,T]$. Denote the mean exposure function by
$\mathbb{E}(X(t))=\mu(t)$. Define the covariance function at any two
timepoints $(s,t)$ as $cov(X(s),X(t))$. Assume that the covariance function
can be expressed by the orthogonal expansion in the $L^{2}$ sense

$cov(X(s),X(t))=\sum_{k=1}^{\infty}\lambda_{k}\phi_{k}(s)\phi_{k}(t)$ (S4)

where
$\lambda_{1}\geqslant\lambda_{2}\geqslant\ldots\geqslant\lambda_{\infty}\geqslant
0$ are called eigenvalues and
$\phi_{1}(\cdot),\phi_{2}(\cdot),\ldots,\phi_{\infty}(\cdot)$ are the
orthonormal eigenfunctions satisfying in sequence

$\phi_{1}=\operatorname*{arg\,max}_{\|\phi\|=1}\left\\{var\left(\int_{0}^{T}(X(t)-\mu(t))\phi(t)dt\right)\right\\}$
$\vdots$
$\phi_{k}=\operatorname*{arg\,max}_{\|\phi\|=1;<\phi,\phi_{j}>=0;j=1,\ldots,k-1}\left\\{var\left(\int_{0}^{T}(X(t)-\mu(t))\phi(t)dt\right)\right\\}$
$\vdots$

where $\|\phi\|:=\left[\int_{0}^{T}\phi(t)^{2}dt\right]^{\frac{1}{2}}$ and
$<\phi,\phi_{j}>:=\int_{0}^{T}\phi(t)\phi_{j}(t)dt$. The $k$-th eigenfunction
is the function maximizing the variance of the projected value of the exposure
curve over this function across individuals and is also orthogonal to the
previous eigenfunctions. We define the principal component associated with the
$k$-th eigenfunction as $\xi_{k}=\int_{0}^{T}(X(t)-\mu(t))\phi_{k}(t)dt$. For
any $k_{1}$-th and $k_{2}$-th principal components, it is easy to know

$\begin{split}cov(\xi_{k_{1}},\xi_{k_{2}})=&cov(\int_{0}^{T}X(t)\phi_{k_{1}}(t)dt,\int_{0}^{T}X(t)\phi_{k_{2}}(t)dt)=cov(\boldsymbol{\phi}_{k_{1}}^{T}\boldsymbol{X}\delta,\boldsymbol{\phi}_{k_{2}}^{T}\boldsymbol{X}\delta)=\delta\boldsymbol{\phi}_{k_{1}}^{T}cov(\boldsymbol{X},\boldsymbol{X})\boldsymbol{\phi}_{k_{2}}\delta\\\
\overset{\left\langle
1\right\rangle}{=}&\delta\boldsymbol{\phi}_{k_{1}}^{T}\left(\sum_{k=1}^{\infty}\lambda_{k}\boldsymbol{\phi}_{k}\boldsymbol{\phi}_{k}^{T}\right)\boldsymbol{\phi}_{k_{2}}\delta=\sum_{k=1}^{\infty}\lambda_{k}(\boldsymbol{\phi}^{T}_{k_{1}}\boldsymbol{\phi}_{k}\delta)(\boldsymbol{\phi}^{T}_{k_{2}}\boldsymbol{\phi}_{k}\delta)\\\
=&\sum_{k=1}^{\infty}\lambda_{k}\int_{0}^{T}\phi_{k_{1}}(t)\phi_{k}(t)dt\int_{0}^{T}\phi_{k_{2}}(t)\phi_{k}(t)dt=0\qquad\textrm{for
}k_{1}\neq k_{2}\end{split}$ (S5)

due to the eigenfunction properties, where
$\boldsymbol{\phi}_{k}=(\phi_{k}(\delta)\,\phi_{k}(2\delta)\,\ldots\,\phi_{k}(T))^{T}$,
$\boldsymbol{X}=(X(\delta)\,X(2\delta)\,\ldots\,X(T))^{T}$ and $\delta$ is a
small time grid for numerical approximation. Equation $\left\langle
1\right\rangle$ is due to the orthogonal expansion (S4).

The variance of principal components in the $k$-th eigenfunction direction is
the $k$-th eigenvalues. That is, $var(\xi_{k})=\lambda_{k}$. According to (S4)
and the orthonormal eigenfunction property, we can know

$\int_{0}^{T}cov(X(s),X(t))\phi_{k}(s)ds=\lambda_{k}\phi_{k}(t)$ (S6)

The curve can be then expressed by the Karhunen-Loéve theorem as
$X(t)-\mu(t)=\sum_{k=1}^{\infty}\xi_{k}\phi_{k}(t)$. Hence, we can
approximately express each individual curve with the first $M$ eigenfunctions
explaining enough variation (e.g., $>95\%$)

$X_{i}(t)\approx\mu(t)+\sum_{m=1}^{M}\xi_{m,i}\phi_{m}(t)$ (S7)

where $\xi_{m,i}$ is the principal component value of the $i$-th individual
w.r.t the $m$-th eigenfunction.

### Text S3: Principal components analysis through conditional expectation

We briefly introduce how to draw FPCA with sparse measured samples. One can
see the original paper yao2005functional proposing PACE in detail. Assume we
have the individual data where the individual measured exposure are
$\\{X_{ij};j=1,\ldots,N_{i}\\}$ where $X_{ij}$ represents the $j$-th
observation at timepoint $t_{ij}$ of the $i$-th individual. PACE makes the
following assumptions:

  * (i)

The measured exposure could contain the measurement error for the accurate
exposure but it is in an additive independent way; that is,
$X_{ij}=X_{i}(t_{ij})+\epsilon_{ij}$ and
$\epsilon_{ij}\overset{\textrm{i.i.d}}{\sim}\mathcal{N}(0,\sigma^{2})$.

  * (ii)

The individual exposure is a random function expressed as the FPCA (via the
Karhunen-Loéve theorem):
$X_{i}(t_{ij})=\mu(t_{ij})+\sum_{k=1}^{\infty}\xi_{k,i}\phi_{k}(t_{ij})$ where
the distribution of the principal components $\\{\xi_{k};k=1,2,\ldots\\}$ are
assumed to be Gaussian.

  * (iii)

The individual measured timepoints $\\{t_{ij};j=1,\ldots,N_{i}\\}$ and the
number of the measured timepoint, $N_{i}$, are independent of the measurement
errors and the principal components.

Therefore, the individual random variables
$\\{\boldsymbol{X}_{i},\boldsymbol{\xi}_{i}\\}$, where
$\boldsymbol{X}_{i}=(X_{i1}\,\cdots\,X_{iN_{i}})^{T}$ and
$\boldsymbol{\xi}_{i}=(\xi_{1,i}\,\cdots\,\xi_{\infty,i})^{T}$, follow the
multivariate normal distribution:

$\left(\begin{array}[]{c}\boldsymbol{\xi}_{i}\\\
\boldsymbol{X}_{i}\end{array}\right)\sim\mathcal{N}(\left(\begin{array}[]{c}\boldsymbol{0}\\\
\boldsymbol{\mu}_{i}\end{array}\right),\left(\begin{array}[]{cc}var(\boldsymbol{\xi}_{i})&\boldsymbol{\Gamma}\\\
\boldsymbol{\Gamma}^{T}&\boldsymbol{\Sigma}_{\boldsymbol{X}_{i}}\end{array}\right))$
(S8)

where
$\boldsymbol{\Gamma}_{k,j}:=cov(\xi_{k,i},X_{i}(t_{ij}))=cov(\xi_{k,i},\xi_{k,i}\phi_{k}(t_{i}j))=\lambda_{k}\phi_{k}(t_{ij})$
due to the FPCA property. $\lambda_{k}=var(\xi_{k})$ is the eigenvalue for the
$k$-th principal component,
$\boldsymbol{\Sigma}_{\boldsymbol{X}_{i}}=var(\boldsymbol{X}_{i})$, and
$\boldsymbol{\mu}_{i}=(\mu(t_{i1})\,\cdots\,\mu(t_{iN_{i}}))^{T}$. The
conditional expectation of the $k$-th principal components given the other
observations for an individual is therefore

$\mathbb{E}(\xi_{k,i}|\boldsymbol{X}_{i})=\lambda_{k}\boldsymbol{\phi}^{T}_{k,i}\boldsymbol{\Sigma}_{\boldsymbol{X}_{i}}^{-1}(\boldsymbol{X}_{i}-\boldsymbol{\mu}_{i})$
(S9)

where
$\boldsymbol{\phi}_{k,i}=(\phi_{k}(t_{i1})\,\cdots\,\phi_{k}(t_{iN_{i}}))^{T}$.
In practice, we first draw the curve and surface smoothing with observed
population data to obtain the estimates $\hat{\mu}(t)$ and
$\widehat{cov}(X(t),X(s))$ (as well as $\hat{\sigma}$ and
$\hat{\boldsymbol{\Sigma}}_{\boldsymbol{X}}$); then derive the estimated
eigenfunctions $\\{\hat{\phi}_{k}(t)\\}$ and eigenvalues
$\\{\hat{\lambda}_{k}\\}$ based on the equation (S6); finally, the conditional
expectation of principal components with the plug-in estimates is

$\hat{\xi}_{k,i}=\hat{\mathbb{E}}(\xi_{k,i}|\boldsymbol{X}_{i})=\hat{\lambda}_{k}\boldsymbol{\hat{\phi}}^{T}_{k,i}\boldsymbol{\hat{\Sigma}}_{\boldsymbol{X}_{i}}^{-1}(\boldsymbol{X}_{i}-\hat{\boldsymbol{\mu}}_{i})$
(S10)

which is used as the $k$-th principal component score of the $i$-th
individual. Similar to before, we can choose the first $K$ principal
components for reasonably good approximation and name $K$ as the number of the
primary principal components. $K$ can be chosen according to the fraction-of-
variance-explained (FVE), cross-validation style, AIC-type criteria, etc.

### Text S4: Two-sample summary-level data

Since we transform the time-varying analysis problem into MVMR problems, one
advantage of this approach is that it allows for a two-sample and summary-
level data setting. This means that the exposure can come from one dataset,
while the outcome can come from another dataset
burgess2015using,davies2018reading. The more flexible data structure
significantly increases the applicability of our methods for time-varying
studies. In most cases, public longitudinal data do not measure both the
exposure and outcome together, as the exposures and outcomes are often defined
after the study protocol. In contrast, the most widely recognized models for
time-varying studies, namely the structural (nested) mean models
robins1989analysis,robins1994correcting,robins2000marginal, typically require
one-sample data for fitting. The corresponding IV models
michael2023instrumental also require complete longitudinal information for the
instrument, exposure, and outcome.

Our MPCMR models can also be implemented with summary-level data. For summary
statistics, we additionally assume that the instrument-exposure model is
(where the principal components are either the original PC or the transformed
ones with transformation matrix, depending on the choice of basis functions):

${\xi}_{k}={\alpha}_{0}+\boldsymbol{\alpha}_{k}^{T}\boldsymbol{Z}+V_{k}\quad
k=1,2,\ldots,K$ (S11)

with $\mathbb{E}(V_{k1}V_{k2}|\boldsymbol{Z})$ does not depend on
$\boldsymbol{Z}$ for any $k_{1}$ and $k_{2}$, which is called the
homoskedasticity assumption. Therefore, we have the instrument-outcome model:
$Y={\theta}_{0}+\boldsymbol{\theta}^{T}\boldsymbol{Z}+U\,,\,k=1,2,\ldots,K$
with $\mathbb{E}(U|\boldsymbol{Z})=0$ and $\mathbb{E}(U^{2}|\boldsymbol{Z})$
does not depend on $\boldsymbol{Z}$. We can construct the moment functions
with summary statistics
$\hat{\boldsymbol{g}}(\boldsymbol{\beta})=\hat{\boldsymbol{\theta}}-\hat{\boldsymbol{\alpha}}\boldsymbol{\beta}$
where the estimated association
$\hat{\boldsymbol{\theta}}\in\mathbb{R}^{J\times 1}$ and
$\hat{\boldsymbol{\alpha}}\in\mathbb{R}^{J\times K}$ can be from either the
univariable regression or multivariable regression. Then the effect function
$\beta(t)$ can be estimated with the similar GMM procedures introduced in the
main text.

### Text S5: MPCMR estimation with IVW

We illustrate the MPCMR estimation with the classical IVW framework. We first
assume the effect shape satisfies the parametric form:
$\beta(t)=\sum_{l=1}^{L}\gamma_{l}b_{l}(t)$ where $\\{\beta_{l}\\}$ are known
basis functions and $L\leq K$ for identification.

Review that the structural equation of MPCMR (5) to form the MVMR respective
to the parameters of interest $\\{\gamma_{l}\\}$

$\begin{split}Y_{i}=&\beta_{0}+\int_{0}^{T}\beta(t)X_{i}(t)dt+g_{Y}(U_{i},\epsilon_{Y,i})\\\
=&\beta_{0}+\sum_{l=1}^{L}\gamma_{l}\int_{0}^{T}b_{l}(t)X_{i}(t)dt+g_{Y}(U_{i},\epsilon_{Y,i})\end{split}$
(S12)

where the MVMR estimates $\hat{\gamma}_{l}$ is based on the IVW regression,
where we assume the SNPs are pruned to be or assumed to be not correlated with
each other,

$\begin{split}\hat{\theta}_{j}=&\sum_{l=1}^{L}\gamma_{l}\frac{cov(G_{j},\int_{0}^{T}b_{l}(t)X(t)dt)}{var(G_{j})}+\epsilon_{j}\qquad\epsilon_{j}\sim\mathcal{N}(0,s.e.(\hat{\theta}_{j})^{2}\tau^{2})\\\
=&\sum_{l=1}^{L}\gamma_{l}\frac{cov(G_{j},\int_{0}^{T}b_{l}(t)\left[\mu(t)+\sum_{k=1}^{K}\xi_{k}\phi_{k}(t)\right]dt)}{var(G_{j})}+\epsilon_{j}\\\
=&\sum_{l=1}^{L}\gamma_{l}\sum_{k=1}^{K}\frac{cov(G_{j},\xi_{k})}{var(G_{j})}\int_{0}^{T}b_{l}(t)\phi_{k}(t)dt+\epsilon_{j}\end{split}$
(S13)

with the plug-in estimates
$\frac{\widehat{cov}(G_{j},\xi_{k})}{\widehat{var}(G_{j})}$ from the
regressing the PCs on the SNPs, and the standard numeric estimates
$\int_{0}^{T}b_{l}(t)\phi_{k}(t)dt\approx\sum_{i}^{I}\phi_{k}(t_{i})b_{l}(t_{i})T/I$
for the uniform grids $\\{t_{i}\\}$. The IVW regression is therefore

$\hat{\boldsymbol{\theta}}=\hat{\boldsymbol{\alpha}}\boldsymbol{B}\boldsymbol{\gamma}+\boldsymbol{\epsilon}\qquad\boldsymbol{\epsilon}\sim\mathcal{N}(\boldsymbol{0},\tau^{2}\boldsymbol{\Sigma})$
(S14)

where
$\boldsymbol{B}=\int_{0}^{T}\boldsymbol{\phi}(t)\boldsymbol{b}^{T}(t)dt\in\mathbb{R}^{K\times
L}$, $\hat{\boldsymbol{\alpha}}\in\mathbb{R}^{J\times K}$ is the instrument-
PCs association estimates and $\boldsymbol{\Sigma}$ is the diagonal matrix
with $\boldsymbol{\Sigma}_{j,j}=s.e.(\hat{\theta}_{j})^{2}$. The estimators
are

$\hat{\boldsymbol{\gamma}}=(\boldsymbol{B}^{T}\hat{\boldsymbol{\alpha}}^{T}\boldsymbol{\Sigma}^{-1}\hat{\boldsymbol{\alpha}}\boldsymbol{B})^{-1}\boldsymbol{B}^{T}\hat{\boldsymbol{\alpha}}^{T}\boldsymbol{\Sigma}^{-1}\hat{\boldsymbol{\theta}}$
(S15)

with the estimator variance matrix
$(\boldsymbol{B}^{T}\hat{\boldsymbol{\alpha}}^{T}\boldsymbol{\Sigma}^{-1}\hat{\boldsymbol{\alpha}}\boldsymbol{B})^{-1}\hat{\tau}^{2}$
and the pointwise functional error
$\boldsymbol{b}^{T}(t)(\boldsymbol{B}^{T}\hat{\boldsymbol{\alpha}}^{T}\boldsymbol{\Sigma}^{-1}\hat{\boldsymbol{\alpha}}\boldsymbol{B})^{-1}\hat{\tau}^{2}\boldsymbol{b}(t)$.
There are typically two common choices for the parametric form of $b_{l}(t)$:

  1. 1.

The polynomial basis function: $b_{l}(t)=t^{l-1}$, $l=1,2,\ldots$

  2. 2.

The eigenfunction system: $b_{l}=\phi_{l}(t)$, $l=1,2,\ldots$

One special kind of parametric assumption used in recent TVMR studies is to
assume the constant-effect form that $\beta(t)=\gamma$ for any timepoint
cao2016mendelian, which leads to a one-parameter problem and greatly
simplifies the estimation model. In our context, such an assumption means
$\boldsymbol{B}=\int_{0}^{T}\boldsymbol{\phi}(t)1dt=\boldsymbol{1}\in\mathbb{R}^{K}$
as the eigenfunction property. The IVW regression becomes the UVMR regression
model:
$\boldsymbol{\theta}_{j}=\left(\sum_{k}^{K}\hat{\alpha}_{j,k}\right)\gamma+\epsilon_{j}$
where $\epsilon_{j}\sim\mathcal{N}(0,s.e.(\theta_{j})^{2})$.

### Text S6: Time-varying IV validity and basis function testing

In Mendelian randomization, some genetic variants may be invalid. This can be
caused by the possible genetic correlation with confounders (hence the
violation of exchangeability), or the genetic pleiotropy phenomenon
davey2003mendelian (hence the violation of the exclusion restriction). One may
wish to assess the IV validity. The IV validity test is also called the over-
identification test hansen1982large, assessing whether all instruments are
estimating the same causal effect parameters ($\\{\beta^{\ast}_{k}\\}$ or
$\\{\gamma_{l}\\}$). One special source making an instrument invalid is the
basis function assumptions for $\beta(t)$. When using the incorrect basis
function in $\boldsymbol{B}$, the instrument for the transformed principal
components will fail to estimate the same parameter and therefore be invalid.
When some instrument is invalid, the over-identification should be rejected.
Note that we cannot test the validity when all instruments used are invalid in
the same way, therefore caution should be given that the IV validity test
assesses only the coherency rather than validity parente2012cautionary.
However, due to the ignorable possibility of all instruments being invalid in
the same way, we can regard the over-identification test as the IV validity
test.

We can test the IV validity using many methods, including Cochran’s Q test
cochran1954combination, Sargan J test sargan1988testing, or other goodness-of-
fit tests. We take the Q test as an example due to its popular usage for
summary-level MR analysis. Assume the null hypothesis $H_{0}$: all instruments
are valid and $\beta(t)=\sum_{l=1}^{L}\gamma_{l}b_{l}(t)$ is the correct shape
form. We assume that each single instrument is uncorrelated with each other.
This independent instrument assumption can be guaranteed in practice by
pruning the SNPs to near independence (e.g. by using the genetics analysis
software PLINK purcell2007plink) or compressing them to principal components
by ordinary PCA batool2022disentangling. When independent instruments cannot
be guaranteed, one can choose the Sargan J test instead (this is also embedded
in our package TVMR). Let the number of instruments and the (untransformed)
PCs be $J$ and $K$ (for identification, we require $J\geq K\geq L$). The
corresponding test statistic is:

$Q=\sum_{j=1}^{J}\frac{(\hat{\theta}_{j}-\hat{\boldsymbol{\alpha}}_{j}^{T}\boldsymbol{B}\hat{\boldsymbol{\gamma}})^{2}}{s.e.(\hat{\theta}_{j})^{2}+\hat{\boldsymbol{\gamma}}^{T}\boldsymbol{B}^{T}\boldsymbol{\Sigma}_{\alpha,j}\boldsymbol{B}\hat{\boldsymbol{\gamma}}-2\hat{\boldsymbol{\gamma}}^{T}\boldsymbol{B}^{T}\boldsymbol{\Gamma}_{j}}$
(S16)

where $\hat{\boldsymbol{\alpha}}_{j}$ is the estimated instrument-PCs
associations for the $j$-th instrument, and $\boldsymbol{\Sigma}_{\alpha,j}$
is the corresponding variance matrix.
$\boldsymbol{\Gamma}_{j}=\text{cov}(\hat{\boldsymbol{\alpha}}_{j},\hat{\theta}_{j})$.
Under the null hypothesis $H_{0}$ that the instruments are valid and the basis
function assumption is correctly specified, the statistic follows a
$\chi^{2}_{J-L}$ distribution asymptotically. The details for calculating each
element in this statistic can be seen in Supplementary Text S9. Particularly,
the quantities $\hat{\boldsymbol{\gamma}}$ can be estimated by either the
maximal profile likelihood procedure zhaoMRRAPS,wang2021causal or the
iterative steps (given in Supplementary Text S9) that are commonly used in
summary-data MR studies bowden2019improving.

Note that one can use the statistic (S16) to test the IV validity only,
regardless of the basis function assumptions. This is because, no matter the
form of $\beta(t)$, the instruments should estimate the same parameters
$\\{\beta^{\ast}_{k}\\}$ under the null hypothesis. This is equivalent to the
classical MVMR model, in which case the Q statistic equals the formula (S16)
by setting $\boldsymbol{B}=\boldsymbol{I}$ or any other transformation matrix
$\boldsymbol{B}\in\mathbb{R}^{K\times K}$ that is full rank.

### Text S7: Time-varying IV strength testing

Another critical problem in Mendelian randomization and instrumental variable
analysis is the weak instrument. The weak instrument refers to the
instrument’s strength being insufficient to support meaningful and reliable
results. A weak instrument can lead to ill-identification problems and make
the inference on the parameters of interest invalid and unreliable
andrews2019weak. When fitting univariable MR models, the strength of the
instrument depends on the strength and uncertainty of the genetic
associations, which can be expressed by its effect size or the F statistics
[ch. 16]hernan2020causal. When fitting multivariable MR models, the concept of
a weak instrument needs to take into account the possible linear dependence
that the genetic association with a single exposure (PC in our context) may be
linearly expressed by the genetic association with the remaining exposures
(PCs). This leads to larger estimation error and, therefore, weak instrument
bias. Such linear dependence can be expressed by either conditional F
statistics sanderson2016weak, sanderson2019examination or the Q statistic
sanderson2021testing. The same test statistic (S16) can be used for testing
the instrument strength.

When using the Q statistic to test the instrument strength for the $k$-th PC
with the remaining PCs, we treat the $k$-th PC as the ’outcome,’ and the
remaining PCs are the ’exposures.’ We then let
$\hat{\theta}_{j}=\hat{\boldsymbol{\alpha}}_{j,k}$ and use the updated
$\hat{\boldsymbol{\alpha}}_{j}\leftarrow\hat{\boldsymbol{\alpha}}_{j,-k}$ with
the corresponding adjusted standard errors. Set
$\boldsymbol{B}=\boldsymbol{I}\in\mathbb{R}^{(K-1)\times(K-1)}$, and the
testing degree of freedom is $J-(K-1)$. Strong evidence for rejection
indicates strong instrument strength.

The IV strength and validity test is provided by our package TVMR.

### Text S8: Related proof of LM statistic

We show that
$\hat{\Delta}_{k}=\hat{\mathbb{E}}_{n}(-\boldsymbol{Z}\xi_{k}(Y-\boldsymbol{\xi}^{T}\boldsymbol{\beta}_{0})\boldsymbol{Z}^{T})$
is the consistent estimator of
$\Delta_{k}:=cov(\sqrt{n}\hat{\boldsymbol{g}}(\boldsymbol{\beta}_{0}),\sqrt{n}\hat{\boldsymbol{G}}_{k})$.

###### Proof.

Without loss of generality, we let all variables be centered; hence

$\begin{split}\Delta_{k}=&cov(\sqrt{n}\hat{\boldsymbol{g}}(\boldsymbol{\beta}_{0}),\sqrt{n}\hat{\boldsymbol{G}}_{k})\\\
=&n\,cov(\hat{\mathbb{E}}_{n}(\boldsymbol{Z}(Y-\boldsymbol{\xi}^{T}\boldsymbol{\beta}_{0})),\hat{\mathbb{E}}_{n}(-\boldsymbol{Z}\xi_{k}))\\\
=&cov(\boldsymbol{Z}(Y-\boldsymbol{\xi}^{T}\boldsymbol{\beta}_{0}),-\boldsymbol{Z}\xi_{k})\quad\textrm{due
to i.i.d samples}\\\
=&\mathbb{E}(-\boldsymbol{Z}\xi_{k}(Y-\boldsymbol{\xi}^{T}\boldsymbol{\beta}_{0})\boldsymbol{Z}^{T})\end{split}$
(S17)

∎

We show that
$cov(\sqrt{n}\hat{D}_{k}(\boldsymbol{\beta}_{0}),\sqrt{n}\hat{\boldsymbol{g}}(\boldsymbol{\beta}_{0}))\overset{p}{\to}\boldsymbol{0}$.

###### Proof.

review that
$\hat{D}_{k}(\boldsymbol{\beta}_{0})=\hat{\boldsymbol{G}}_{k}-\hat{\Delta}_{k}(\boldsymbol{\beta}_{0})^{T}\hat{\Omega}(\boldsymbol{\beta}_{0})^{-1}\hat{\boldsymbol{g}}(\boldsymbol{\beta}_{0})$;
hence

$\begin{split}cov(\sqrt{n}\hat{D}_{k}(\boldsymbol{\beta}_{0}),\sqrt{n}\hat{\boldsymbol{g}}(\boldsymbol{\beta}_{0}))=&cov(\sqrt{n}\hat{\boldsymbol{G}}_{k},\sqrt{n}\hat{\boldsymbol{g}}(\boldsymbol{\beta}_{0}))-cov(\sqrt{n}\hat{\Delta}_{k}(\boldsymbol{\beta}_{0})^{T}\hat{\Omega}(\boldsymbol{\beta}_{0})^{-1}\hat{\boldsymbol{g}}(\boldsymbol{\beta}_{0}),\sqrt{n}\hat{\boldsymbol{g}}(\boldsymbol{\beta}_{0}))\\\
\overset{p}{\to}&cov(\sqrt{n}\hat{\boldsymbol{G}}_{k},\sqrt{n}\hat{\boldsymbol{g}}(\boldsymbol{\beta}_{0}))-{\Delta}_{k}(\boldsymbol{\beta}_{0})^{T}{\Omega}(\boldsymbol{\beta}_{0})^{-1}cov(\sqrt{n}\hat{\boldsymbol{g}}(\boldsymbol{\beta}_{0}),\sqrt{n}\hat{\boldsymbol{g}}(\boldsymbol{\beta}_{0}))\\\
=&\Delta_{k}(\boldsymbol{\beta}_{0})^{T}-{\Delta}_{k}(\boldsymbol{\beta}_{0})^{T}{\Omega}(\boldsymbol{\beta}_{0})^{-1}{\Omega}(\boldsymbol{\beta}_{0})\quad\textrm{by
definition}\\\ =&\boldsymbol{0}\end{split}$ (S18)

∎

### Text S9: Q statistic testing details

We provide the testing inference details regarding the parametric MPCMR
models. For UVMR, a similar inference procedure can be found in the papers
zhaoMRRAPS,bowden2019improving; for MVMR, see the papers
wang2021causal,sanderson2021testing. Assume the null hypothesis $H_{0}$ is:
$\beta(t)=\sum_{l=1}^{L}\gamma_{l}b_{l}(t)$ for fixed basis functions
$\\{b_{l}(t)\\}$ and parameters $\\{\gamma\\}$ to be estimated. We assume that
the SNPs are pruned to be uncorrelated to each other; therefore, the
covariance of the genetic associations from different SNPs to any phenotypes
(PCs, outcome, whether same or different) is neglectably zero. We can build
the joint likelihood for the estimated associations with the phenotypes for
each SNP

$\left(\begin{array}[]{c}\hat{\theta}_{j}\\\
\hat{\boldsymbol{\alpha}}_{j}\end{array}\right)\sim\mathcal{N}(\left(\begin{array}[]{c}\boldsymbol{\alpha}_{j}^{T}\boldsymbol{B}\boldsymbol{\gamma}\\\
\boldsymbol{\alpha}_{j}\end{array}\right),\left(\begin{array}[]{cc}s.e.(\hat{\theta}_{j})^{2}&\boldsymbol{\Gamma}_{j}^{T}\\\
\boldsymbol{\Gamma}_{j}&\boldsymbol{\Sigma}_{\alpha,j}\\\
\end{array}\right))\qquad j=1,2,\ldots,J$ (S19)

Where $\boldsymbol{\Sigma}_{\alpha,j}$ is the variance matrix for
$\hat{\boldsymbol{\alpha}}_{j}$. We can approximate that
$\boldsymbol{\Sigma}_{\alpha,j}$ is a diagonal matrix with the $k$-th element
equal to $s.e.(\hat{\alpha}_{j,k})^{2}$.
$\boldsymbol{\Gamma}_{j}:=cov(\hat{\boldsymbol{\alpha}}_{j},\hat{\theta}_{j})$.
When the exposure and the outcome come from two samples,
$\boldsymbol{\Gamma}_{j}=\boldsymbol{0}$; otherwise $\boldsymbol{\Gamma}_{j}$
is estimated by the approximation under the overlapping-samples setting

$\boldsymbol{\Gamma}_{j,k}=cov(\hat{\alpha}_{j,k},\hat{\theta}_{j})=\frac{n_{s}}{n_{1}n_{2}}\frac{1}{var(G_{j})}cov(v^{PC}_{k},v^{Y})$
(S20)

with the relevant plug-in estimates, where $v^{PC}_{k}$ and $v^{Y}$ are the
residuals of the regression on the genetic variants from the $k$-th PC and the
outcome, respectively, under the common dataset; $n_{1}$ and $n_{2}$ are the
sample size of the PC data and the outcome data, $n_{s}$ is the shared sample
size.

###### Proof.

Assume the genetic variants are uncorrelated with each other. Consider the
regression for the two phenotypes $\xi_{k}$ and $Y$ with the respective
(possibly overlapping) samples:

$\xi_{k,i}=\alpha_{0}+\boldsymbol{\alpha}^{T}\boldsymbol{G}_{i}+v_{1,i}\qquad
i=1,\ldots,n_{1}$ (S21)
$Y_{i}=\theta_{0}+\boldsymbol{\theta}^{T}\boldsymbol{G}_{i}+v_{2,i}\qquad
i=1,\ldots,n_{2}$ (S22)

where $v_{1}\overset{\textrm{def}}{=}v^{PC}_{k}$ and
$v_{2}\overset{\textrm{def}}{=}v^{Y}$;
$\alpha_{j}=cov(G_{j},\xi_{k})/var(G_{j})$ and
$\theta_{j}=cov(G_{j},Y)/var(G_{j})$. When $\boldsymbol{G}$ are binary (like
in most MR cases), it is easy to know $\mathbb{E}(G^{p}v_{1})=0$ and
$\mathbb{E}(G^{p}v_{2})=0$ for any order $p=0,1,\ldots$. Therefore, the
genetic association estimators for the $j$-th genetic variant are

$\hat{\alpha}_{j,k}=\frac{\widehat{cov}(G_{j},\xi_{k})}{\widehat{var}(G_{k})}=\frac{(\overline{G_{j}\xi_{k}})_{n_{1}}-(\overline{G_{j}})_{n_{1}}(\overline{\xi_{k}})_{n_{1}}}{\widehat{var}(G_{k})}$
(S23)
$\hat{\theta}_{j}=\frac{\widehat{cov}(G_{j},Y)}{\widehat{var}(G_{k})}=\frac{(\overline{G_{j}Y})_{n_{2}}-(\overline{G_{j}})_{n_{2}}(\overline{Y})_{n_{2}}}{\widehat{var}(G_{k})}$
(S24)

The approximation has the following steps

$\begin{split}cov(\hat{\alpha}_{j,k},\hat{\theta}_{j})\overset{\left\langle
1\right\rangle}{=}&\frac{1}{var(G_{j})^{2}}cov[(\overline{G_{j}\xi_{k}})_{n_{1}}-(\overline{G_{j}})_{n_{1}}(\overline{\xi_{k}})_{n_{1}},(\overline{G_{j}Y})_{n_{1}}-(\overline{G_{j}})_{n_{1}}(\overline{Y})_{n_{1}}]\\\
=&\frac{1}{var(G_{j})^{2}}cov[\frac{1}{n_{1}}\sum_{i=1}^{n_{1}}G_{j,i}(\xi_{k,i}^{C}+(\overline{\xi_{k}})_{n_{1}})-(\overline{G_{j}})_{n_{1}}(\overline{\xi_{k}})_{n_{1}},\frac{1}{n_{2}}\sum_{i=1}^{n_{1}}G_{j,i}(Y_{i}^{C}+(\overline{Y})_{n_{2}})-(\overline{G_{j}})_{n_{2}}(\overline{Y})_{n_{2}}]\\\
\overset{\left\langle
2\right\rangle}{=}&\frac{1}{var(G_{j})^{2}}cov[\frac{1}{n_{1}}\sum_{i=1}^{n_{1}}G_{j,i}\xi_{k,i}^{C},\frac{1}{n_{2}}\sum_{i=1}^{n_{2}}G_{j,i}Y_{i}^{C}]\\\
\overset{\left\langle
3\right\rangle}{=}&\frac{n_{s}}{n_{1}n_{2}}\frac{1}{var(G_{j})^{2}}cov(G_{j}\xi_{k}^{C},G_{j}Y^{C})\\\
\overset{\left\langle
4\right\rangle}{=}&\frac{n_{s}}{n_{1}n_{2}}\frac{1}{var(G_{j})^{2}}cov(\alpha_{j,k}G^{2}_{j}+G_{j}v_{1},\theta_{j}G^{2}_{j}+G_{j}v_{2})\\\
=&\frac{n_{s}}{n_{1}n_{2}}\frac{1}{var(G_{j})^{2}}\left(\alpha_{j,k}\theta_{j}var(G_{j}^{2})+cov(G_{j}v_{1},G_{j}v_{2})\right)\\\
\overset{\left\langle
5\right\rangle}{\approx}&\frac{n_{s}}{n_{1}n_{2}}\frac{1}{var(G_{j})^{2}}cov(G_{j}v_{1},G_{j}v_{2})\\\
\overset{\left\langle
6\right\rangle}{=}&\frac{n_{s}}{n_{1}n_{2}}\frac{1}{var(G_{j})}cov(v_{1},v_{2})\end{split}$
(S25)

Where the centred phenotypes $\xi_{k}^{C}$ and $Y^{C}$ in Equation
$\left\langle 2\right\rangle$ corresponds to the equations
$\xi_{k}-\mathbb{E}(\xi_{k})=\boldsymbol{\alpha}^{T}\boldsymbol{G}+v_{1}$ and
$Y-\mathbb{E}(Y)=\boldsymbol{\theta}^{T}\boldsymbol{G}+v_{2}$, and $n_{s}$ in
Equation $\left\langle 3\right\rangle$ is the shared sample size. In Equation
$\left\langle 4\right\rangle$, we assume
$\mathbb{E}(G^{2}_{j}G_{j_{\ast}}v_{1})=0$ and
$\mathbb{E}(G^{2}_{j}G_{j_{\ast}}v_{2})=0$ for any $j_{\ast}$. We approximate
$\alpha_{j,k}\theta_{j}\approx 0$ in Equation $\left\langle 5\right\rangle$,
as the genetic associations in Mendelian randomization are usually small. In
equation $\left\langle 6\right\rangle$, we assume
$\mathbb{E}(G_{j}^{2}v_{1}v_{2})=\mathbb{E}(G_{j}^{2})\mathbb{E}(v_{1}v_{2})$,
which equals to $var(G_{j})cov(v_{1},v_{2})$. ∎

Hence, the estimator for $\boldsymbol{\Gamma}_{j,k}$ is

$\hat{\boldsymbol{\Gamma}}_{j,k}=\widehat{cov}(\hat{\alpha}_{j,k},\hat{\theta}_{j})=\frac{n_{s}}{n_{1}n_{2}}\frac{1}{\widehat{var}(G_{j})}\widehat{cov}(v^{PC}_{k},v^{Y})\approx\frac{1}{n_{1}n_{2}}\frac{1}{\widehat{var}(G_{j})}\sum_{i=1}^{n_{s}}(\xi_{k,i}-\hat{\xi}_{k,i})(Y_{i}-\hat{Y}_{i})\qquad
k=1,\ldots,K$ (S26)

where $\widehat{cov}(v^{PC}_{k},v^{Y})$ is estimated by the share sample part,
and we can approximate that
$\frac{1}{n_{s}}\sum_{i=1}^{n_{s}}(\xi_{k,i}-\hat{\xi}_{k,i})\sum_{i=1}^{n_{s}}(Y_{i}-\hat{Y}_{i})\approx
0$. The statistic is therefore

$Q=\sum_{j=1}^{J}\frac{(\hat{\theta}_{j}-\hat{\boldsymbol{\alpha}}_{j}^{T}\boldsymbol{B}\hat{\boldsymbol{\gamma}})^{2}}{s.e.(\hat{\theta}_{j})^{2}+\hat{\boldsymbol{\gamma}}^{T}\boldsymbol{B}^{T}\boldsymbol{\Sigma}_{\alpha,j}\boldsymbol{B}\hat{\boldsymbol{\gamma}}-2\hat{\boldsymbol{\gamma}}^{T}\boldsymbol{B}^{T}\boldsymbol{\Gamma}_{j}}$
(S27)

Under the null that the parametric model is correctly specified, the statistic
follows the chi-squared distribution $\chi^{2}_{J-L}$. The estimator of
$\hat{\boldsymbol{\gamma}}$ based on the statistic is regarded as the robust
estimator. The iterative steps for deriving the robust estimator and Q
statistic is

  1. (1)

Obtain the initial estimates $\hat{\boldsymbol{\gamma}}_{(0)}$ based on
Equation (S15), which is to minimize

$Q(\boldsymbol{\gamma};\boldsymbol{\gamma}_{\ast})=\sum_{j=1}^{J}\frac{(\hat{\theta}_{j}-\hat{\boldsymbol{\alpha}}_{j}^{T}\boldsymbol{B}\boldsymbol{\gamma})^{2}}{s.e.(\hat{\theta}_{j})^{2}+\boldsymbol{\gamma}_{\ast}^{T}\boldsymbol{B}^{T}\boldsymbol{\Sigma}_{\alpha,j}\boldsymbol{B}\boldsymbol{\gamma}_{\ast}-2\boldsymbol{\gamma}_{\ast}^{T}\boldsymbol{B}^{T}\boldsymbol{\Gamma}}$
(S28)

with $\boldsymbol{\gamma}_{\ast}=\boldsymbol{0}$ (this is equivalent to the
standard MVMR fitting).

  2. (2)

Minimize $Q(\boldsymbol{\gamma};\hat{\boldsymbol{\gamma}}_{(0)})$ to obtain
$\hat{\boldsymbol{\gamma}}_{(1)}$

  3. (3)

Repeat the above step and update the iterative estimates
$\hat{\boldsymbol{\gamma}}_{(n)}$ until the estimates are stable.

The corresponding $Q$ statistic value with the plug-in MPCMR robust estimates
is then used to evaluate the parametric shape assumptions or other testing
objectives (like IV strength and validity as introduced in the previous
supplementary texts).

plain All_TVMR