# Robust inference for the unification of confidence intervals in meta-
analysis

Wei Liang, Haicheng Huang, Xiamen University, China,
Hongsheng Dai∗, Newcastle University, Newcastle upon Tyne, U.K.
<EMAIL_ADDRESS>
and Yinghui Wei, University of Plymouth, U.K.

###### Abstract

Traditional meta-analysis assumes that the effect sizes estimated in
individual studies follow a Gaussian distribution. However, this
distributional assumption is not always satisfied in practice, leading to
potentially biased results. In the situation when the number of studies,
denoted as $K$, is large, the cumulative Gaussian approximation errors from
each study could make the final estimation unreliable. In the situation when
$K$ is small, it is not realistic to assume the random-effect follows Gaussian
distribution. In this paper, we present a novel empirical likelihood method
for combining confidence intervals under the meta-analysis framework. This
method is free of the Gaussian assumption in effect size estimates from
individual studies and from the random-effects. We establish the large-sample
properties of the non-parametric estimator, and introduce a criterion
governing the relationship between the number of studies, $K$, and the sample
size of each study, $n_{i}$. Our methodology supersedes conventional meta-
analysis techniques in both theoretical robustness and computational
efficiency. We assess the performance of our proposed methods using simulation
studies, and apply our proposed methods to two examples. confidence Interval;
Empirical Likelihood; Meta-analysis; Random-effect Model.

00footnotetext: To whom correspondence should be addressed.

##  1 Introduction

Meta-analysis is a statistical method to combine results from multiple
studies, facilitating the inference-making process concerning specific
parameters of interests, such as odds ratio. Parameters estimated from
individual studies can be combined using either a fixed-effect or a random-
effects meta-analysis model (Higgins and Thompson, 2009; Deeks _and others_ ,
2022). A fixed-effect model assumes that a single parameter value is common to
all studies, and a random-effects model assumes that parameters of the
underlying studies follow some distributions.

Let the random variable $Y_{i}$ $(i=1,2,\cdots,K)$ denote an effect size
estimate from study $i$. Here $Y_{i}$ could be the mean difference between
treatment groups, log-odds ratios, log-hazard ratios, or other effect size
estimates from study $i$. The conventional random-effects meta-analysis model
is defined as

$Y_{i}=\theta_{i}+\epsilon_{i},\;\;\;\;\;\theta_{i}=\theta+\xi_{i}$ (1)

where $\theta_{i}$ is the true effect size of study $i$, $\theta$ is the
common overall effect size, $\epsilon_{i}$ is the within-study random error
and $\xi_{i}$ is a random variable reflecting study-specific deviation from
the common overall effect. In the literature, it is often assumed that
$\epsilon_{i}$ and $\xi_{i}$ are independent Gaussian variables, with
$\epsilon_{i}\sim{\cal N}(0,\sigma_{i}^{2})$ and $\xi_{i}\sim{\cal
N}(0,\tau^{2})$, where $\sigma_{i}^{2}$ is the within-study error variance and
depends on the sample size $n_{i}$ of study $i$. For instance,
$\sigma_{i}^{2}$ is usually in the order of $n_{i}^{-1}$. The between-study
variance $\tau^{2}$ reflects the heterogeneity in the effect sizes across
studies. Under the conventional random-effects meta-analysis model, the
marginal distribution of $Y_{i}$ is a Gaussian distribution with mean $\theta$
and variance $\sigma_{i}^{2}+\tau^{2}$. The overall effect $\theta$ can be
estimated by the maximum likelihood method (Brockwell and Gordon, 2001;
Normand, 1999); and the between-study variance $\tau^{2}$ can also be
estimated by the (profile) maximum likelihood method (Brockwell and Gordon,
2001; Hardy and Thompson, 1996; Partlett and Riley, 2017; Langan _and others_
, 2019), DerSimonian and Laird’s method of moments (DerSimonian and Laird,
1986) or quantile approximation (Brockwell and Gordon, 2007). Equation (1)
becomes a fixed-effect meta-analysis model when $\xi_{i}=0$
($\mathbb{E}\,\xi_{i}=0$ and $\tau^{2}=0$). Within the existing literature,
both models have leaned upon a Gaussian assumption. Nonetheless, practical
reality often deviates from this parametric assumption, as demonstrated by
(Lee and Thompson, 2008; Brockwell and Gordon, 2001; Sutton and Higgins,
2008).

One important objective in meta-analysis is to estimate a confidence interval
for the overall effect size. The confidence interval can be constructed using
the asymptotic distribution for the estimator $\hat{\theta}$ , derived from a
maximum likelihood approach (Brockwell and Gordon, 2001) or a profile
likelihood approach (Brockwell and Gordon, 2001; Hardy and Thompson, 1996).
Both approaches relied on the assumption that the estimator $\hat{\theta}$
approximately follows a Gaussian distribution, which however may not be true
if $n_{i}$ and/or $K$ is small. Constructing a confidence interval for
$\theta$ based on alternative parametric distributions was proposed, including
$t$-distribution, other parametric distributions which allow for potentially
skewed and heavy tails in random-effects distributions (Sidik and Jonkman,
2002; Lee and Thompson, 2008; Baker and Jason, 2008; Ozturk and Balakrishnan,
2020), and quantile approximation (Brockwell and Gordon, 2007; Jackson and
Bowden, 2009). Confidence distributions (CDs), which can be viewed as
“distribution estimators” (Schweder and Hjort, 2002), is a technique which can
be used to construct point estimators, confidence intervals and $p$-values.
The fusion of CD functions from individual studies have been explored for
constructing a confidence interval for the overall effect size (Singh _and
others_ , 2005; Xie _and others_ , 2011). Nevertheless, the CD approach still
relies on Gaussian assumptions for the integration of CD functions from
individual studies into a unified CD function. For example, Singh _and others_
(2005) and Xie _and others_ (2011) used Gaussian distribution to construct CD
functions and Nagashima _and others_ (2019) proposed a parametric bootstrap
algorithm based on Gaussian and $t$-distributions. We refer the reader to
Veroniki _and others_ (2018) for a review of methods (all depend on parametric
assumptions) to calculate a confidence interval for the estimated overall
effect size from a random-effects meta-analysis.

Different from the existing meta-analysis methods, we consider a more
challenging situation where only the $1-\alpha$ confidence interval for the
effect size of study $i$, $[L_{i},U_{i}]$, is observed, and individual study
effect size estimates $Y_{i}$ and within-study error variance estimate
$\hat{\sigma}_{i}^{2}$ are not observed. Further, we seek to relax the
Gaussian assumption. The aim of this paper is to propose a novel non-
parametric meta-analysis method based on empirical likelihood (EL) to make
robust inferences on the confidence interval of the overall effect size
$\theta$. We relax the Gaussian distributional assumption for $Y_{i}$ and
$\theta_{i}$. In other words, we consider

$Y_{i}=\theta_{i}+\epsilon_{i},\;\;\epsilon_{i}\sim
F_{i}(\cdot);\qquad\theta_{i}=\theta+\xi_{i},\;\;\xi_{i}\sim G(\cdot),$ (2)

where the distributions $F_{i},G$ are not subject to any parametric
assumption. The proposed non-parametric method allows $K\rightarrow\infty$,
which is more reasonable than many existing literature where $K$ was treated
as a constant. This innovation empowers the non-parametric meta-analysis to
effectively manage a substantial number of studies.

The rest of the paper is organised as follows. The proposed non-parametric
meta-analysis method is presented in Section 2 and 3. We evaluate the
performance of the proposed method using simulation studies in Section 4. Two
real example applications are presented in Section 5. The paper concludes with
a discussion and theoretical proofs for the large sample theory are provided
in Web Appendix A of the Supplementary Material.

##  2 Model-based meta-analysis approach

For the parametric Gaussian models presented in Section 1, Singh _and others_
(2005) and Xie _and others_ (2011) developed a general framework for combining
results from individual studies in a meta-analysis. They used a Confidence
Distribution (CD) function $H(\bm{X},\theta)$ of the data $\bm{X}$ and
parameter $\theta$, from which inference on confidence intervals, point
estimates or $p$-values can be derived. Their method amalgamates the
individual study’s CD functions $H_{i}(\bm{X}_{i},\theta)$ into a unified CD
function $H_{c}(\theta)$, through a monotonic function
$g_{c}(H_{1}(\bm{X}_{1},\theta),\cdots,H_{K}(\bm{X}_{K},\theta))$. Under the
Gaussian assumption in the random-effects meta-analysis model, and assuming
the estimated effect size $Y_{i}$ and the within-study error variance
$\sigma_{i}^{2}$ for studiy $i$ are available, a typical CD function for
$\theta$ is

$H_{i}(\theta)=\Phi\left(\frac{\theta-Y_{i}}{s_{i}}\right),$ (3)

where $\Phi(\cdot)$ is the cumulative Gaussian distribution function,
$s_{i}^{2}={\sigma}_{i}^{2}+\hat{\tau}^{2}$ in the random-effects model
($s_{i}^{2}={\sigma}_{i}^{2}$ in the fixed-effect model), and $\hat{\tau}^{2}$
is the estimated between-study variance.

Singh _and others_ (2005) and Xie _and others_ (2011) summarized methodologies
used in the literature in group decision analysis, where linear combination of
individual study result is commonly used as the final outcome. They introduced
the following linear combination function

$\displaystyle
g_{c}(u_{1},\cdots,u_{K})=\sum_{k=1}^{K}w_{k}\,\Phi^{-1}(u_{k}),$ (4)

where the weight $w_{k}=1/s_{k}$. Then the combination of CDs for $\theta$ is
defined as

$\displaystyle H_{c}(\theta)$
$\displaystyle=\Phi\left(\frac{g_{c}(H_{1}(\theta),\cdots,H_{K}(\theta))}{\sqrt{\sum_{i=1}^{K}w_{i}^{2}}}\right)=\Phi\left(\frac{\sum_{k=1}^{K}w_{k}\,\frac{\theta-
Y_{k}}{s_{k}}}{\sqrt{\sum_{i=1}^{K}w_{i}^{2}}}\right).$

Inference on $\theta$ can then be obtained using the combined CD function
$H_{c}(\theta)$. For instance, the $(1-\beta)$ confidence interval for
$\theta$ can be constructed as

$[L_{C},\,U_{C}]=[H_{c}^{-1}(\beta/2),\,H_{c}^{-1}(1-\beta/2)].$ (5)

The above approach involves two assumptions: 1\. The model error terms in
model (1) in Section 1 follow a Gaussian distribution; 2. The linear
combination in (4) follows a zero-mean Gaussian distribution. If each of the
sample size $n_{i}$ is large enough, the model error term in model (1) can be
assumed to be Gaussian, i.e. the above assumption 1 is reasonable, but it
fails if $n_{i}$ is small. However, this is a strong assumption. If the first
assumption is violated, the second assumption will also be false. If the
Gaussian assumption is not held true for some individual studies, in the
presence of a large number of studies, i.e. large $K$ , the accumulated errors
by combining all different studies might not be bounded. For instance, under
model (2) without the Gaussian assumptions, if we consider $\theta_{i}\sim
N(\theta=0,1)$ but
$y_{ij}|\theta_{i}\sim\chi^{2}(4)-4+\theta_{i},Y_{i}=n_{i}^{-1}\sum_{j=1}^{n_{i}}y_{ij}$,
then we will not be able to use the combined statistic
$Z=K^{-\frac{1}{2}}\sum_{i=1}^{K}(\theta-Y_{i})/s_{i}$ to construct the meta-
analysis confidence interval, since $Z$ will not converge to $0$ when
$K\rightarrow\infty$. This can be seen from Figure 1. Furthermore, the
relationship between the number of studies $K$ and sample size $n_{i}$ in
meta-analysis has not been extensively discussed. We will propose a non-
parametric method and relax this Gaussian assumption in the following section
and discuss the relationship between $K$ and $n_{i}$.

Figure 1: The empirical mean of $Z$ with $\theta_{i}\sim
N(0,1),\,y_{ij}|\theta_{i}\sim\chi^{2}(4)-4+\theta_{i}$ and
$Y_{i}=n_{i}^{-1}\sum_{j=1}^{n_{i}}y_{ij}$. We consider the setting with large
$K=100,500,1000,5000,10000,50000$ and small $n_{i}=20,40,60$ for all $i$.

##  3 Non-parametric meta-analysis based on empirical likelihood

In this section, we propose a robust inference procedure for conducting meta-
analysis that does not make distributional assumptions on effect size
estimates for individual studies. Hence, this new approach is more reliable
even if some studies have a small sample size $n_{i}$. We first present a
generalized theorem and necessary assumptions.

Let $\bm{W}_{i}(\theta)$ be an $r$-dimensional function of $\theta$ for study
$i$. If $\bm{W}_{i}(\theta)$ is asymptotically unbiased to $0$ at the true
parameter value $\theta_{0}$, we can construct the estimating equation
$\frac{1}{K}\sum_{i=1}^{K}\bm{W}_{i}(\theta)=0$ for the unknown parameter
$\theta$. Following standard conventions in the literature, individual studies
are independent of each other. Hence, we assume that the following assumptions
hold.

######  Assumption 1.

Asymptotic Unbiaseness. For any $i$, at the true parameter $\theta_{0}$,

${\mathbb{E}}\,\bm{W}_{i}(\theta_{0})=O(n_{i}^{-1/2}),$

which means that each component of ${\mathbb{E}}\,\bm{W}_{i}$ is in the order
of $n_{i}^{-1/2}$.

######  Assumption 2.

Independence. $\bm{W}_{i}(\theta),i=1,\cdots,K$ are independent.

We then consider the Empirical Likelihood (EL) approach and define the EL
function

$\mathcal{L}(\theta)=\sup\left\\{\prod_{i=1}^{K}p_{i}:\,\sum_{i=1}^{K}p_{i}\,\bm{W}_{i}(\theta)=0,\,\sum_{i=1}^{K}p_{i}=1\right\\}.$

At the true parameter $\theta_{0}$, the corresponding EL ratio statistic is

$\displaystyle\mathcal{R}(\theta_{0})=\frac{\mathcal{L}(\theta_{0})}{\sup_{\theta}\mathcal{L}(\theta)}.$
(6)

We further need the following assumptions to discuss the large-sample
properties of the EL ratio statistic.

######  Assumption 3.

Finite Variance. For any $i$, at the true parameter $\theta_{0}$,
$\bm{W}_{i}(\theta_{0})$ has a full rank finite covariance matrix
$\textbf{M}_{i}=\text{var }(\bm{W}_{i}(\theta_{0}))$, and
$\textbf{M}_{i},i=1,\cdots,K$ are uniformly bounded by M, which means for any
vector $\textbf{e}\in\mathbb{R}^{r}$,
$\sup_{i}\textbf{e}^{{\mathrm{\scriptscriptstyle
T}}}\,\textbf{M}_{i}\,\textbf{e}\leq\textbf{e}^{{\mathrm{\scriptscriptstyle
T}}}\,\textbf{M}\,\textbf{e}$.

Note that Assumption 3 always holds in practice, since typical classical
statistical estimation can always be formalised as a sequence of
asymptotically unbiased estimating equations, such as the likelihood inference
and generalised estimating equation approaches. The finite variance of
$\textbf{M}_{i}$ is guaranteed when we use standardised estimating equations.

The last assumption is about the relationship between the number of studies
$K$ and the individual study sample size $n_{i}$, which is necessary when we
discuss the large sample theories of meta-analysis under a non-parametric
approach. To illustrate this assumption, we introduce the following notation
and assumption. Let $n^{(K)}$ be the minimum number of individual sample
sizes, that is $n^{(K)}=\min\\{n_{i},i=1,\cdots,K\\}$.

######  Assumption 4.

Large Sample Order Condition. $K/n^{(K)}\rightarrow 0$, as
$K\rightarrow\infty$ or $n^{(K)}\rightarrow\infty$.

Assumption 4 includes the scenario in the literature, where $K$ is treated as
a fixed number and sample sizes $n_{i}$ may goes to infinity. First, the above
assumption allows $K$ to be a fixed number and in this case, as all $n_{i}$
goes to infinity, the above assumption is always true. In addition, the above
assumption is more general since it allows allows $K$, the number of studies,
to go to infinity. Assumption 4 implies that for large values of $K$, we only
need $K$ to be in a lower order than all $n_{i}$, the sample sizes of
individual studies. It gives a sufficient condition required for that meta-
analyses can be conducted with a substantial multitude of studies. Recently,
Cai _and others_ (2023) discussed the case of large $K$ but finite $n_{i}$.

Under these assumptions, we have the main result as presented in the following
theorem, with the proofs given in Web Appendix A of the Supplementary
Material.

######  Theorem 5.

Assume that $\mathcal{R}(\theta)$ is an EL ratio statistic defined in (6).
Under Assumptions 1-4, as $K\rightarrow\infty$, we have
$-2\log\mathcal{R}(\theta_{0})\rightarrow\chi^{2}(r)$ in distribution.

For a fixed confidence level $\beta$, we can construct the $(1-\beta)$-level
confidence interval for $\theta$ as

$\displaystyle\\{\theta:-2\log\mathcal{R}(\theta)\leq\chi^{2}_{1-\beta}(r)\\}.$
(7)

Next, we will apply this theorem to two different meta-analysis models. The
key point in applying the EL inference is to construct an estimating equation
to define the EL ratio statistic.

###  3.1 Random-effects meta-analysis

Consider the random-effects meta-analysis model (2). Assume that we have
collected $1-\alpha$ confidence intervals $[L_{i},\,U_{i}]$ for the parameter
$\theta_{i}$ from $K$ studies. The sample size for each study is denoted as
$n_{i}$, $i=1,\,2,\,\cdots,\,K$. Based on the confidence intervals
$[L_{i},\,U_{i}]$, we define

$W_{i}(\theta)=\frac{U_{i}+L_{i}-2\theta}{U_{i}-L_{i}},$ (8)

which can be written as
$W_{i}(\theta)=\frac{\overline{G}_{i}(\theta)-\underline{G}_{i}(\theta)}{\overline{G}_{i}(\theta)+\underline{G}_{i}(\theta)}$,
with $\overline{G}_{i}(\theta)=U_{i}-\theta$ and
$\underline{G}_{i}(\theta)=\theta-L_{i}$. Then

$\displaystyle{\mathbb{E}}(\underline{G}_{i}(\theta_{0}))$ $\displaystyle=$
$\displaystyle{\mathbb{E}}\,(\theta_{0}-Y_{i}+Y_{i}-L_{i})={\mathbb{E}}\,({\mathbb{E}}(\theta_{0}-Y_{i})|\theta_{i})+{\mathbb{E}}(Y_{i}-L_{i})$
$\displaystyle=$
$\displaystyle{\mathbb{E}}(\theta_{0}-\theta_{i})+{\mathbb{E}}(Y_{i}-L_{i})={\mathbb{E}}(Y_{i}-L_{i}),$

and ${\mathbb{E}}(\overline{G}_{i}(\theta_{0}))={\mathbb{E}}(U_{i}-Y_{i})$.
Because of the large sample property of $Y_{i}$, it is generally true that

$\displaystyle{\mathbb{E}}(Y_{i}-L_{i})\,=\,O(n_{i}^{-1/2}),\qquad{\mathbb{E}}(U_{i}-Y_{i})\,=\,O(n_{i}^{-1/2}).$

On the other hand, when using empirical likelihood method to construct
confidence intervals, we have
${\mathbb{E}}(Y_{i}-L_{i})-{\mathbb{E}}(U_{i}-Y_{i})\,\leq c\cdot n_{i}^{-1}$,
for some constant $c$. This holds as well when confidence intervals are
constructed via symmetric distributions, such as t-distribution or Gaussian
distribution, since ${\mathbb{E}}(Y_{i}-L_{i})-{\mathbb{E}}(U_{i}-Y_{i})=0$.
Therefore, we here focus on the confidence intervals such that,

$\displaystyle{\mathbb{E}}(Y_{i}-L_{i})-{\mathbb{E}}(U_{i}-Y_{i})\leq c\cdot
n_{i}^{-1}.$

Now we have

${\mathbb{E}}(\overline{G}_{i}(\theta_{0}))-{\mathbb{E}}(\underline{G}_{i}(\theta_{0}))\,=\,O(n_{i}^{-1}),\qquad{\mathbb{E}}(\overline{G}_{i}(\theta_{0}))+{\mathbb{E}}(\underline{G}_{i}(\theta_{0}))\,=\,O(n_{i}^{-1/2}),$
(9)

and thus Assumption 1 is satisfied, i.e.,

${\mathbb{E}}\,W_{i}(\theta_{0})=O(n_{i}^{-1/2}).$ (10)

Now we define

$\displaystyle\mathcal{L}(\theta)$
$\displaystyle=\sup\left\\{\prod_{i=1}^{K}p_{i}:\,\sum_{i=1}^{K}p_{i}\left(\frac{L_{i}+U_{i}-2\theta}{U_{i}-L_{i}}\right)=0,\,\sum_{i=1}^{K}p_{i}=1\right\\},$

then the corresponding EL ratio statistic is
$\mathcal{R}(\theta_{0})=\mathcal{L}(\theta_{0})/\sup_{\theta}\mathcal{L}(\theta).$
According to Theorem 5, we have
$-2\log\mathcal{R}(\theta_{0})\rightarrow\chi^{2}(1)$ in distribution.

###  3.2 Fixed-effect meta-analysis

We now consider a fixed-effect model, i.e. for all $i=1,\ldots,K$,
$\theta_{i}=\theta$ in equation (2). Since this is a special case of the
random-effects model, we can still use the method presented in the previous
subsection to make inferences. In the fixed-effect model, we can also use the
confidence level to construct the estimating equation. First, we define

$\displaystyle V_{i}(\theta)=I\\{L_{i}\leq\theta\leq U_{i}\\}-(1-\alpha).$
(11)

Most of the time, the confidence interval $[L_{i},\,U_{i}]$ is constructed
based on an asymptotic distribution rather than an accurate distribution,
hence, the expectation of (11) is approximately equal to $0$. For example,
when the confidence interval $[L_{i},U_{i}]$ was constructed according to the
central limit theorem, we have

$\displaystyle{\mathbb{E}}\,V_{i}(\theta_{0})={\mathbb{E}}\,(I\\{L_{i}\leq\theta_{0}\leq
U_{i}\\}-(1-\alpha))=O(n_{i}^{-1/2}).$

However, in the random-effects meta-analysis model the above
$V_{i}(\theta_{0})$ is not asymptotically unbiased to $0$. This is
demonstrated by the following result under model (1) with the Gaussian
assumption for $\epsilon_{i}$ and $\xi_{i}$,

$\displaystyle{\mathbb{E}}\,V_{i}(\theta_{0})$
$\displaystyle={\mathbb{E}}\,\left(I\left\\{Y_{i}-z_{1-\frac{\alpha}{2}}\cdot\hat{\sigma}_{i}\leq\theta_{0}\leq
Y_{i}+z_{1-\frac{\alpha}{2}}\cdot\hat{\sigma}_{i}\right\\}\right)-(1-\alpha)$
$\displaystyle={\mathbb{E}}\,\left(\Phi\left(z_{1-\frac{\alpha}{2}}-\frac{\theta_{i}-\theta_{0}}{\hat{\sigma}_{i}}\right)-\Phi\left(-z_{1-\frac{\alpha}{2}}-\frac{\theta_{i}-\theta_{0}}{\hat{\sigma}_{i}}\right)\right)-(1-\alpha).$

This expectation does not converge to $0$, and it cannot be used to construct
the asymptotically unbiased estimating equation in the random-effects model.
Therefore, we consider the following three ways to define EL ratio statistic
in the fixed-effect meta-analysis model. Denote $\bm{V}_{i}(\theta)$ as the
estimating equation for study $i$, and $r$ as the dimension of
$\bm{V}_{i}(\theta)$.

  * •

Case I: the estimating equation based on the confidence level, i.e. the
following function $V_{i}(\theta)=I\\{L_{i}\leq\theta\leq
U_{i}\\}-(1-\alpha)$.

  * •

Case II: the estimating equation based on symmetry (this is the same as the
random-effects model), i.e. using the following function

$V_{i}(\theta)=W_{i}(\theta)=z_{1-\frac{\alpha}{2}}\,\frac{L_{i}+U_{i}-2\theta}{U_{i}-L_{i}}.$

  * •

Case III: the estimating equation based on confidence level and symmetry. In
this case, $r=2$, and
$V_{i}(\theta)=(V_{i1}(\theta),\,\,V_{i2}(\theta))^{{\mathrm{\scriptscriptstyle
T}}}$, where

$V_{i1}(\theta)=I\\{L_{i}\leq\theta\leq U_{i}\\}-(1-\alpha),\quad
V_{i2}(\theta)=z_{1-\frac{\alpha}{2}}\,\frac{L_{i}+U_{i}-2\theta}{U_{i}-L_{i}}.$

By choosing an appropriate $V_{i}(\theta)$, we can define three different
empirical likelihood statistics, which allow us to make inferences on the
confidence interval in meta-analysis,

$\displaystyle\mathcal{L}_{f1}(\theta)$
$\displaystyle=\sup\left\\{\prod_{i=1}^{K}p_{i}:\,\sum_{i=1}^{K}p_{i}I\\{L_{i}\leq\theta\leq
U_{i}\\}=1-\alpha,\,\sum_{i=1}^{K}p_{i}=1\right\\},$
$\displaystyle\mathcal{L}_{f2}(\theta)$
$\displaystyle=\sup\left\\{\prod_{i=1}^{K}p_{i}:\,\sum_{i=1}^{K}p_{i}\left(\frac{L_{i}+U_{i}-2\theta}{U_{i}-L_{i}}\right)=0,\,\sum_{i=1}^{K}p_{i}=1\right\\},$
(12) $\displaystyle\mathcal{L}_{f3}(\theta)$
$\displaystyle=\sup\left\\{\prod_{i=1}^{K}p_{i}:\,\sum_{i=1}^{K}p_{i}I\\{L_{i}\leq\theta\leq
U_{i}\\}=1-\alpha,\,\sum_{i=1}^{K}p_{i}\left(\frac{L_{i}+U_{i}-2\theta}{U_{i}-L_{i}}\right)=0,\,\sum_{i=1}^{K}p_{i}=1\right\\}.$

Then the corresponding empirical likelihood ratio statistics are

$\displaystyle\mathcal{R}_{fk}(\theta_{0})=\frac{\mathcal{L}_{fk}(\theta_{0})}{\sup_{\theta}\mathcal{L}_{fk}(\theta)},\qquad
k=1,2,3$

respectively. Following Theorem 5, we have
$-2\log\mathcal{R}_{f1}(\theta_{0})\rightarrow\chi^{2}(1)$,
$-2\log\mathcal{R}_{f2}(\theta_{0})\rightarrow\chi^{2}(1)$ and
$-2\log\mathcal{R}_{f3}(\theta_{0})\rightarrow\chi^{2}(2)$, in distribution.
Then for a fixed confidence level $\beta$, we can construct the
$(1-\beta)$-level confidence interval for $\theta$ based on (7).

##  4 Simulation Studies

To demonstrate the performance of the proposed non-parametric method, we
conduct a comprehensive comparison against two established methods: the
existing CD method and the conventional meta-analysis method. To be specific,
DerSimonian & Laird method is used to estimate $\tau^{2}$ in CD method, and
‘rma’ in the ‘metafor’ R package where we use REML to estimate $\tau^{2}$.

We consider four different scenarios in this section, which corresponds to
model (2) but have different distributions for the zero-mean error terms,

$\displaystyle y_{ij}$ $\displaystyle\sim
F(\cdot\mid\theta_{i},\sigma^{2}),\;\;Y_{i}=n_{i}^{-1}\sum_{j=1}^{n_{i}}y_{ij},$
$\displaystyle\theta_{i}$ $\displaystyle=\theta+\xi_{i},\;\;\;\xi_{i}\sim
G(\cdot\mid\tau^{2}),$ (13)

where $\mathbb{E}(y_{ij})=\theta_{i}$, i.e. $\theta_{i}$ is the expected
effect size of study $i$, $Y_{i}$ is the estimated effect size of study $i$,
and $\theta$ is the overall effect size. The parameter $\sigma^{2}$ is the
within-study variance. This model setting is the same as (2), where
var$(Y_{i})$ is in the order of $\sigma_{i}^{2}=\sigma^{2}/n_{i}$. For all
scenarios, we set $\theta=0$, and $\sigma^{2}=1$, and we consider different
values of $K$, $20,60$ and $100$.

We set the four simulation scenarios as follows, with different distributions
$F$ and $G$.

  * •

Scenario 1 is the Gaussian distribution case, where $F={\cal
N}(\theta_{i},\sigma^{2})$ and $G={\cal N}(0,\tau^{2})$. This is equivalent to
model (1).

  * •

Scenario 2 sets $F={\cal N}(\theta_{i},\sigma^{2})$ and $G$ to log-normal
distribution with mean $0$ and variance $\tau^{2}$. The individual study
effect $\theta_{i}$ is not Gaussian. This is common in practice. For example,
a study’s effect size may have a positive and heavy tail. We choose
$\xi_{i}=c_{2}\cdot\text{log-normal}(1,1)+d_{2}$, where $c_{2}$ and $d_{2}$
are constants such that the distribution of $\xi_{i}$ has mean $0$ and
variance as a pre-specified value $\tau^{2}$.

  * •

Scenario 3 sets $F$ to be non-Gaussian, but $G$ to be Gaussian ${\cal
N}(0,\tau^{2})$. When an individual study has a small sample size, the
individual studies’ effect size estimate $Y_{i}$ does not follow a Gaussian
distribution. In our simulation, we first generate $\theta_{i}\sim{\cal
N}(0,\tau^{2}),\,i=1,\,2,\cdots,\,K$. Given the sample size $n_{i}$ and
$\theta_{i}$, draw $y_{ij}\mid\theta_{i}\sim a_{3,i}\cdot\chi^{2}(4)+b_{3,i}$
such that
${\mathbb{E}}(y_{ij}\mid\theta_{i})=\theta_{i},\text{var}(y_{ij}\mid\theta_{i})=\sigma^{2}$.

  * •

Scenario 4 is the most general case where we set both $F$ and $G$ to be non-
Gaussian. We choose $G$ to be the distribution of
$c_{4}\cdot\chi^{2}(4)+d_{4}$ to generate $\theta_{i}$, and we obtain random
sample $Y_{i}=n_{i}^{-1}\sum_{j=1}^{n_{i}}y_{ij}$ with
$y_{ij}\mid\theta_{i}\sim a_{4,i}\cdot\text{log-normal}(1,1)+b_{4,i}$, where
$a_{4,i}$ and $b_{4,i}$ are such that
${\mathbb{E}}(y_{ij}\mid\theta_{i})=\theta_{i},$ and
$\text{var}(y_{ij}\mid\theta_{i})=\sigma^{2}$.

For each scenario, we obtain the estimated effect size and the estimated
within-study variance by $Y_{i}=n_{i}^{-1}\sum_{j=1}^{n_{i}}y_{ij}$ and
$\hat{\sigma}_{i}^{2}=\sum_{j=1}^{n_{i}}(y_{ij}-Y_{i})^{2}/(n_{i}(n_{i}-1))$,
respectively. Then the $1-\alpha=95\%$ two-side confidence interval
$[L_{i},U_{i}]$ of $\theta$ is generated by either

$\displaystyle L_{i}=Y_{i}-z_{1-\alpha/2}\cdot\hat{\sigma}_{i},\quad
U_{i}=Y_{i}+z_{1-\alpha/2}\cdot\hat{\sigma}_{i},\qquad n_{i}\geq 30,$ (14)

for large $n_{i}$, or

$\displaystyle L_{i}=Y_{i}-t_{1-\alpha/2}(n_{i}-1)\cdot\hat{\sigma}_{i},\quad
U_{i}=Y_{i}+t_{1-\alpha/2}(n_{i}-1)\cdot\hat{\sigma}_{i},\qquad n_{i}<30,$
(15)

for small $n_{i}$, where $t_{\beta}(n)$ is the $\beta$-quantile for the
$t$-distribution with degree of freedom $n$.

We consider the settings of different values of $\tau^{2}$, ranging from $0$
to $100$. When $\tau^{2}=0$, it becomes a fixed-effect meta-analysis model.
Based on the observations $[L_{i},\,U_{i}]$, we compare $1-\beta=95\%$
confidence intervals which are constructed via the conventional meta-analysis
method, the CD method and the non-parametric EL method.

###  4.1 Scenarios 1, 2 and 3

We consider the conventional settings of random-effects meta-analysis model
$(\tau^{2}\neq 0)$ and fixed-effect meta-analysis model $(\tau^{2}=0)$. Figure
2 shows the results coverage probability against the ratio of between-study
variance and within-study variance for Scenarios 1 to 3, based on different
values of $K$, 20, 60, and 100, $\sigma^{2}=1$, and different values of
$\tau^{2}$, $0,0.001,0.01,0.1,1,10$ and $100$. For each study, the number of
observations $n_{i}$ is drawn randomly from the uniform distribution
$0.2\cdot{\cal U}[100,500]$.

Figure 2: The coverage probabilities of $1-\beta=95\%$ confidence intervals
which are constructed via the conventional meta-analysis method, CD method and
EL method respectively under the setting of different $K$ and $\tau^{2}$ in
Scenario 1. In each setting the $n_{i}$ is sampled from $0.2\cdot{\cal
U}[100,500]$.

The simulation results for Scenario 1 are shown in the first row of plots in
Figure 2. For all three methods the coverage probabilities of the confidence
interval for the overall effect are close to the nominal level $0.95$, under
all different values of $\tau^{2}$. Under the Gaussian assumptions, the non-
parametric EL method is as good as the conventional meta-analysis and the CD
method.

In Scenario 2, the individual study’s effect size estimate is Gaussian but the
random-effect is a log-normal distribution. From the second row of plots in
Figure 2, we can see that all three methods, the conventional meta-analysis,
the CD approach and the non-parametric EL method give underestimated coverage
probabilities, for large values of $\tau^{2}$ (with values $1,10,100$) and
small values of $K$. When $K=20$, the larger $\tau^{2}$ is, the worse
performance for the conventional meta-analysis and the CD method. This is
because the Gaussian assumption for the individual study’s effect size is
invalid and the error is mainly from the non-Gaussian random-effects. When
$\tau^{2}$ becomes larger, the error of using a Gaussian distribution to
approximate a log-normal distribution becomes larger, especially when $K$ is
small. The performance of all three methods improves as $K$ increases. For
example. when $K=100$, all three methods provide reasonably well coverage
probability for overall effect estimate.

The proposed non-parametric EL method exhibits a notable and significant
advantage over the other methods in Scenario 3, where the study effect size is
a $\chi^{2}$ distribution but the random-effect is a Gaussian distribution.
The results are presented in the third row of plots in Figure 2. The coverage
probabilities of confidence intervals derived from both the conventional and
CD methods exhibit a rapid decline as the number of studies ($K$) increases,
especially when the true $\tau^{2}$ value is small. This decline implies a
notable departure from the anticipated large sample properties for these
methods.

In cases where $\tau^{2}$ is small, it is understood that the primary source
of error originates from the within-study approximation. Given that the true
distribution for each study in Scenario 3 follows a $\chi^{2}$ distribution,
the within-study error, under the Gaussian assumption, becomes considerably
inflated. Consequently, the conventional and CD methods manifest
ineffectiveness. However, as the $\tau^{2}$ value increases, errors in the
overall effect estimation stem mainly from the random effects. This scenario
proves beneficial for the conventional and CD methods, as the Gaussian
approximation aligns more closely with the true distribution of $\theta_{i},$
resulting in smaller between-study approximation errors. In contrast, the non-
parametric EL method is more robust and perform better for both small and
large values of $\tau^{2}$. Just as with the conventional and CD methods,
there is a potential for convergence challenges within the non-parametric EL
method. Notably, as $K$ increases, the convergence probability fails to align
with the true nominal value of $0.95$. This deviation can be attributed to the
insufficiency of the sample size ($n_{i}$) within each study, as governed by
Assumption 4. Further exploration of the effect of augmenting the sample size
($n_{i}$) will be undertaken in the subsequent subsection, focusing on
Scenario 4.

###  4.2 Scenario 4 and results based on different $n_{i}$

We now study the performance of different methods under Scenario 4, that both
the individual study effect size and random-effects are non-Gaussian. The top-
row plots within Figure 3 present results for relatively small $n_{i}$ values,
drawn from $0.2\cdot\mathcal{U}[100,500].$ It’s evident that neither the
conventional meta-analysis method nor the CD method performs well under small
$\tau^{2}$ values. Although the non-parametric EL method exhibits an improved
coverage probability, it still falls considerably short of the target nominal
value of $0.95$. This performance gap can be attributed to the relatively
modest $n_{i}$ values presented in this scenario. Significantly larger values
of $n_{i}$ are needed to guarantee the convergence of the log-likelihood ratio
statistic. The bottom row of plots in Figure 3 shows the performance with
large values of $n_{i}$, which are random values from ${\cal U}[4K,5K]$. We
can see that the coverage probabilities based on the conventional meta-
analysis method and the CD method converge very slowly to the true nominal
value $0.95$, even if we have a large number of $n_{i}$ and a large number of
$K$. For instance, when $\tau^{2}=0,K=100$ and $n_{i}$ is drawn from ${\cal
U}[400,500]$, both the conventional method and the CD method give coverage
probabilities about $0.4$. Remarkably, the non-parametric EL method exhibits
significantly improved coverage probabilities, outperforming the conventional
and CD methods. Under the same conditions, the non-parametric EL method
approximates a coverage probability of about $0.8,$ converging much faster to
the true nominal value.

Figure 3: The coverage probabilities of $1-\beta=95\%$ confidence intervals
which are constructed via conventional meta-analysis method, CD method and EL
method respectively under the setting of different $K$ and $\tau^{2}$ in
Scenario 4. In each setting the $n_{i}$ is sampled from
$0.2\cdot\mathcal{U}[100,500]$ (top panel) or sampled from ${\cal U}[4K,5K]$
(bottom panel).

###  4.3 The relationship between $K$ and $n_{i}$ in the Scenario 4

Further, we study the relationship between $K$ and $n_{i}$ in different
methods under Scenario 4 with fixed $\tau^{2}=0.01$. The simulation results
are show in Figure 4. When $K$ is fixed, the coverage probability of each
method will converge to the true nominal value as $n_{i}$ increases, but our
EL method performs better and robustly especially in the small $n_{i}$. When
$n_{i}$ is fixed, the coverage probability of each method will not converge as
$K$ increases since Assumption 4 is not satisfied.

Figure 4: The coverage probabilities of $1-\beta=95\%$ confidence intervals
which are constructed via conventional meta-analysis method, CD method and EL
method respectively under the setting of different $K$ and $n_{i}$ in Scenario
4. In each setting the $n_{i}$ is the same for each $i$ and $\tau^{2}=0.01$.

###  4.4 Convergence of the non-parametric EL estimator

Assumption 4 was not well addressed in the meta-analysis literature, but it
was the key condition to guarantee the consistency of the meta-analysis
estimate. Here we further use a numerical study to demonstrate the convergence
of the distribution of the log-likelihood function at the true parameter,
under different settings of $K$ and $n_{i}$, i.e. numerical justification of
Theorem 5. For brevity, here we only consider the convergence of different
methods mentioned in the fixed-effect model. We use QQ plots to empirically
examine the accuracy of the $\chi^{2}(r)$ approximation to the distribution of
the log-likelihood function at the true parameter,
$-2\log\mathcal{R}(\theta_{0})$.

Taking $F(\cdot\,\mid\,\theta_{i})=\chi^{2}(4),\,\theta_{i}=\theta_{0}=4$ for
the fixed-effect meta-analysis model, and $n_{1}=n_{2}=\cdots=n_{K}=n$, we
compare the sample quantiles obtained from $1000$ replicated values of
$-2\log\mathcal{R}(\theta_{0})$ under different methods. The theoretical
asymptotic distributions of $-2\log\mathcal{R}(\theta_{0})$ under EL1, EL2 and
EL3 are $\chi^{2}(r)$ distributions, with $r=1$ for EL1 and EL2, and $r=2$ for
EL3. The sample quantiles are plotted in Figure 5. The proximity of sample
quantiles to the 45-degree straight line in QQ plots indicates the degree of
resemblance between the distribution of these samples and the $\chi^{2}(r)$
distribution. From these QQ plots we can conclude that, when $K/n=o(1)$,
except for extreme upper quantiles, the distributions of
$-2\log\mathcal{R}(\theta_{0})$ under different methods are reasonably well
approximated by a $\chi^{2}(r)$ distribution. A specific instance can be
observed in the top-right panel of Figure 5, where $n=800$ and $K=200$. Here,
the alignment between sample quantiles and theoretical quantiles is strikingly
close. When $K/n$ increases, especially when $K=800$ and $n=200$, the
convergence of $-2\log\mathcal{R}(\theta_{0})$ to $\chi^{2}(r)$ diminishes. It
is shown by the left-bottom plot of Figure 5, that
$-2\log\mathcal{R}(\theta_{0})$ does not converge to $\chi^{2}(r)$, since
Assumption 4 may be violated. When $n=K$ (the top-left plot and the bottom-
right plot), $-2\log\mathcal{R}(\theta_{0})$ seems still to converge to
$\chi^{2}(r)$, but its performance is not as good as the result with $n=800$
and $K=200$.

Figure 5: A QQ plot with sample quantiles obtained from $1000$ replicated
values of $-2\log\mathcal{R}(\theta_{0})$ under different methods (EL1, EL2
and EL3) and the corresponding the quantiles of the $\chi^{2}(r)$ distribution
($r=1$ for EL1 and EL2, $r=2$ for EL3) under different settings:
$(n,K)=(200,200),(200,800),(800,200),(800,800)$. The 45-degree line is
included in each plot.

##  5 Application

We apply our proposed methods to two real examples. Example 1 has a large
number of studies ($K=74$), whilst Example 2 has a small number of studies
($K=6$). The assumption for $K$ to be bounded in the existing methodologies
may not be valid for a large value $K$ such as in Example 1, since the
accumulated errors from each Gaussian approximation within each study may be
too large. For a small number of $K$ such as in Example 2, the Gaussian
assumption for $\theta_{i}$ in the conventional meta-analysis methods may not
be valid.

###  5.1 Example 1: Mortality in patients with non-COVID illnesses during and
before pandemic

During the COVID-19 pandemic, healthcare providers struggled to manage and
deliver quality health care services. Quantifying the effect of COVID-19
pandemic on patients with non-COVID illnesses and on healthcare systems
compared to pre-pandemic times can help to understand the impact and
disruption of pandemic on health care. A systematic review and meta-analysis
of mortality in patients with non-COVID illness during and before the COVID-19
pandemic was conducted (Lau _and others_ , 2022), which presented a meta-
analysis (74 observational studies with 491,862 patients), and showed an
increased mortality rate during COVID-19 pandemic compared with a pre-
pandenmic period for patients with non-COVID illnesses. They studied the log-
risk ratio (RR), where

$\text{risk ratio}=\frac{\text{risk during the pandemic}}{\text{risk before
the pandemic}},$

where risk is defined as the ratio of the number of death and the number of
patients; and the log-odds ratio (OR), where

$\text{odds ratio}=\frac{\text{odds during the pandemic}}{\text{odds before
the pandemic}},$

where the odds is the ratio of the number of death and the number of alive.

Table 1 shows the confidence intervals for the estimated overall effect. We
also used the the CD function method and our proposed non-parametric empirical
likelihood method. The pandemic is associated with an increased risk of
mortality, with consistent results across the conventional meta-analysis, CD
and EL methods as well as fixed-effect and random-effects meta-analysis. Our
EL approach resulted in larger CIs ($\mbox{RR:}\;1.206\;(1.104,1.469)$,
$\mbox{OR:}\;1.537,\;(1.293,2.070)$), which is expected as our method is non-
parametric.

Due to the heterogeneity of different studies, a random-effects model is more
realistic. Table 1 shows that, when the CD method is used, the risk ratio is
1.384 (95%CI: 1.278 to 1.500) from the random-effects model, and is attenuated
to 1.075 (95%CI: 1.068 to 1.082) when we used a fixed-effect model. However,
when using our EL method, the risk ratios are closer to each other between the
random-effects and fixed-effect meta-analyses. In particular, the non-
parametric approach with empirical likelihood ${\cal L}_{f2}$ in (3.2) can be
used in both random-effects model and fixed-effect model. The results are
similar despite under different model assumptions. We observed similar
patterns for OR in this example.

Table 1: Example 1. The combined confidence interval under each method based on the confidence interval of log-risk ratio or log-odds ratio in each study reported in Lau _and others_ (2022). Combined confidence intervals were estimated by the conventional meta-analysis, CD method, EL methods in random-effects model and fixed-effect model. The conventional meta-analysis was fitted by using ‘rma’ in the ‘metafor’ R, and for the conventional random-effect meta-analysis we used REML method to estimate the between-study variance $\tau^{2}$. | Fixed-effect meta-analysis | Random-effects meta-analysis
---|---|---
| Method | Est. | 95% CI | Method | Est. | 95% CI
Risk Ratio | Conventional | 1.075 | $(1.068,1.082)$ | Conventional | 1.455 | $(1.270,1.668)$
CD | 1.075 | $(1.068,1.082)$ | CD | 1.384 | $(1.278,1.500)$
EL1 | 1.166 | $(1.086,1.443)$ | EL | 1.206 | $(1.104,1.469)$
EL2 | 1.206 | $(1.104,1.469)$ |  |  |
EL3 | 1.175 | $(1.117,1.447)$ |  |  |
Odds Ratio | Conventional | 1.318 | $(1.284,1.353)$ | Conventional | 1.647 | $(1.407,1.929)$
CD | 1.318 | $(1.284,1.353)$ | CD | 1.650 | $(1.401,1.944)$
EL1 | 1.397 | $(1.385,1.547)$ | EL | 1.537 | $(1.293,2.070)$
EL2 | 1.537 | $(1.293,2.070)$ |  |  |
EL3 | 1.423 | $(1.385,1.577)$ |  |  |

###  5.2 Example 2: Stillbirths in pregnancy people with and without COVID-19
vaccination

During the pandemic, COVID-19 vaccination was developed and rolled out at an
unprecedented scale and speed addressing the global public health emergency.
However, COVID-19 vaccination hesitancy has been identified as an issue.
Apprehension surrounding pregnancy was identified as one contributing source
of vaccination hesitancy (Bhattacharya _and others_ , 2022; Golder _and
others_ , 2023; Razai _and others_ , 2021). Evidence of the effectiveness of
COVID-19 vaccination in pregnant women can help to improve the uptake of the
vaccination in this vulnerable group. Stillbirth is an important perinatal
outcome, which is the death or loss of a baby before or during delivery. The
impact of COVID-19 mRNA vaccination on stillbirth was evaluated in a
systematic review (Prasad _and others_ , 2022) and a meta-analysis of six
studies (66,067 vaccinated v.s. 424,624 unvaccinated pregnant people).

Table 2 showed that results are inconsistent across methods. From the results
obtained from the CD method, COVID-19 vaccination is associated with lower
risk of stillbirths. However, in other methods (conventional random-effects
meta-analysis, EL random-effects meta-analysis and EL fixed-effect meta-
analysis), we found no evidence for the association of COVID-19 vaccination
with stillbirths. In summary, the CD method implies that the COVID-19
vaccination is associated with a reduced risk of still-birth. However, the
results from meta-analysis using the non-parametric EL method and the
conventional meta-analysis method cannot confirm this.

Table 2: Example 2. Stillbirths in pregnancy people with and without COVID-19 vaccination. confidence intervals were estimated by using the conventional random-effects and fixed-effect meta-analysis, the CD method and different non-parametric empirical likelihood methods (EL), based on the confidence interval of log-risk ratio or log-odds ratio in each study reported in Prasad _and others_ (2022). | Fixed-effect meta-analysis | Random-effects meta-analysis
---|---|---
| Method | Est. | 95% CI | Method | Est. | 95% CI
Risk Ratio | Conventional | 0.859 | $(0.739,0.997)$ | Conventional | 0.859 | $(0.737,1.001)$
CD | 0.859 | $(0.739,0.997)$ | CD | 0.859 | $(0.739,0.997)$
EL1 | 0.644 | $(0.644,1.163)$ | EL | 0.867 | $(0.762,1.030)$
EL2 | 0.867 | $(0.762,1.030)$ |  |  |
EL3 | 0.854 | $(0.736,1.080)$ |  |  |
Odds Ratio | Conventional | 0.858 | $(0.738,0.997)$ | Conventional | 0.859 | $(0.736,1.001)$
CD | 0.858 | $(0.738,0.997)$ | CD | 0.858 | $(0.738,0.997)$
EL1 | 0.643 | $(0.643,1.163)$ | EL | 0.867 | $(0.761,1.030)$
EL2 | 0.867 | $(0.761,1.030)$ |  |  |
EL3 | 0.856 | $(0.735,1.080)$ |  |  |

##  6 Discussion

This paper introduces a novel non-parametric method for conducting meta-
analysis of confidence intervals. The approach is rooted in empirical
likelihood, applicable to both fixed-effect and random-effects meta-analysis
models. Notably, the non-parametric EL method departs from the conventional
Gaussian assumption often relied upon, while still yielding large sample
theoretical outcomes. Our contribution provides the solutions in addressing
the requirements regarding the number of studies ($K$) and the sample size
($n_{i}$) within each study, to ensure the consistency and asymptotic
normality of the overall effect estimate. This theoretical result fills a
crucial gap that has not been addressed in the current landscape of meta-
analysis literature.

The distinguish feature of the non-parametric EL approach is in its reliance
on an asymptotically unbiased estimating equation for the unknown parameter,
derived from the confidence intervals obtained from individual studies. As
long as the function $\bm{W}_{i}(\theta)$ within the estimating equation
$\bm{W}_{i}(\theta)=\bm{0}$ remains asymptotically unbiased at a rate
$\mathbb{E}\bm{W}_{i}(\theta_{0})=O(n_{i}^{-1/2})$ (as per Assumption 1), the
non-parametric EL approach guarantees a consistent estimate. This assumption,
being quite mild, holds true for a wide array of existing statistical
estimators, such as the maximum likelihood estimator derived from score
functions, the estimating equations in Cox partial likelihood for survival
analysis, and the generalized estimating equation approach.

In our illustrative examples, our non-parametric EL method produced wider
95%CI for the overall effect. In example 2, the confidence intervals estimated
from the EL method are consistent across methods, in contrast, confidence
intervals estimated by CD methods or the conventional meta-analysis are not
consistent between the fixed-effect random-effect models. In particular, it
reveals distinctive insights that diverge from existing meta-analysis
techniques. In Example 2, the non-parametric EL method leads to different
conclusion - there is insufficient evidence that COVID-19 vaccination is
associated with the risk of stillbirths. This exemplifies the effectiveness of
the proposed approach in addressing the underlying challenge in meta-analysis,
particularly when individual study estimates violate the Gaussian assumptions.

## References

  * Baker and Jason (2008) Baker, Rose and Jason, Dan. (2008).  A new approach to outliers in meta-analysis.  Health Care Management Science 11, 121–131.
  * Bhattacharya _and others_ (2022) Bhattacharya, Oashe, Siddiquea1, Bodrun Naher, Shetty, Aishwarya, Afroz, Afsana and Billah, Baki. (2022).  Covid-19 vaccine hesitancy among pregnant women: a systematic review and meta-analysis.  BMJ Open.
  * Brockwell and Gordon (2001) Brockwell, S.E. and Gordon, I.R. (2001).  A comparison of statistical methods for meta-analysis.  Statistics in Medicine 20(16), 825–840.
  * Brockwell and Gordon (2007) Brockwell, Sarah E. and Gordon, Ian R. (2007).  A simple method for inference on an overall effect in meta-analysis.  Statistics in Medicine 26(25), 4531–4543.
  * Cai _and others_ (2023) Cai, Chencheng, Chen, Rong and Xie, Min-ge. (2023).  Individualized group learning.  Journal of the American Statistical Association 118(541), 622–638.
  * Deeks _and others_ (2022) Deeks, Jonathan J., Riley, Richard D. and Higgins, Julian P.T. (2022).  Combining Results Using Meta-Analysis, Chapter 9. John Wiley & Sons, Ltd, pp. 159–184.
  * DerSimonian and Laird (1986) DerSimonian, R. and Laird, N. (1986).  Meta-analysis in clinical trials.  Controlled Clinical Trials 7(3), 177–188.
  * Golder _and others_ (2023) Golder, S., McRobbie-Johnson, A. C. E., Klein, A., Polite, F. G. and G., Gonzalez Hernandez. (2023).  Social media and covid-19 vaccination hesitancy during pregnancy: a mixed methods analysis.  Scandinavian Journal of Statistics 29, 309–332.
  * Hardy and Thompson (1996) Hardy, R.J. and Thompson, S.R. (1996).  A likehood approach to meta-analysis with random effects.  Statistics in Medicine 15(6), 619–629.
  * Higgins and Thompson (2009) Higgins, J.P.T. and Thompson, S.G. (2009).  A re-evaluation of random-effects meta-analysis.  Journal of Royal Statistical Society Series A 172(1), 137–159.
  * Jackson and Bowden (2009) Jackson, Dan and Bowden, Jack. (2009).  A re-evaluation of the ‘quantile approximation method’ for random effects meta-analysis.  Statistics in Medicine 28(2), 338–348.
  * Langan _and others_ (2019) Langan, Dean, Higgins, Julian P.T., Jackson, Dan, Bowden, Jack, Veroniki, Areti Angeliki, Kontopantelis, Evangelos, Viechtbauer, Wolfgang and Simmonds, Mark. (2019).  A comparison of heterogeneity variance estimators in simulated random-effects meta-analyses.  Research Synthesis Methods 10(1), 83–98.
  * Lau _and others_ (2022) Lau, Vincent Issac, Dhanoa, Sumeet, Cheema, Harleen, Lewis, Kimberley, Geeraert, Patrick, Lu, David, Merrick, Benjamin, Vander Leek, Aaron, Sebastianski, Meghan, Kula, Brittany, Chaudhuri, Dipayan, Agarwal, Arnav, Niven, Daniel J., Fiest, Kirsten M., Stelfox, Henry T., Zuege, Danny J., Rewa, Oleksa G. _and others_. (2022, June).  Non-COVID outcomes associated with the coronavirus disease-2019 (COVID-19) pandemic effects study (COPES): A systematic review and meta-analysis.  PLOS ONE 17(6), e0269871.
  * Lee and Thompson (2008) Lee, Katherine J. and Thompson, Simon G. (2008).  Flexible parametric models for random-effects distributions.  Statistics in Medicine 27(3), 418–434.
  * Nagashima _and others_ (2019) Nagashima, Kengo, Noma, Hisashi and Furukawa, Toshi A. (2019).  Prediction intervals for random-effects meta-analysis: A confidence distribution approach.  Statistical Methods in Medical Research 28(6), 1689–1702.
  * Normand (1999) Normand, S.L. (1999).  Meta-analysis: formulating, evaluating, combining, and reporting.  Statistics in Medicine 18(3), 321–359.
  * Ozturk and Balakrishnan (2020) Ozturk, Omer and Balakrishnan, Narayanaswamy. (2020).  Meta-analysis of quantile intervals from different studies with an application to a pulmonary tuberculosis data.  Statistics in Medicine 39(29), 4519–4537.
  * Partlett and Riley (2017) Partlett, Christopher and Riley, Richard D. (2017).  Random effects meta-analysis: Coverage performance of 95% confidence and prediction intervals following reml estimation.  Statistics in Medicine 36(2), 301–317.
  * Prasad _and others_ (2022) Prasad, Smriti, Kalafat, Erkan, Blakeway, Helena, Townsend, Rosemary, O’Brien, Pat, Morris, Edward, Draycott, Tim, Thangaratinam, Shakila, Le Doare, Kirsty, Ladhani, Shamez, Von Dadelszen, Peter, Magee, Laura A., Heath, Paul _and others_. (2022, May).  Systematic review and meta-analysis of the effectiveness and perinatal outcomes of COVID-19 vaccination in pregnancy.  Nature Communications 13(1), 2414.
  * Razai _and others_ (2021) Razai, Mohammad S, Chaudhry, Umar A R, Doerholt, Katja, Bauld, Linda and Majeed, Azeem. (2021).  Covid-19 vaccination hesitancy.  BMJ 373.
  * Schweder and Hjort (2002) Schweder, T. and Hjort, N. L. (2002).  Confidence and likelihood.  Scandinavian Journal of Statistics 29, 309–332.
  * Sidik and Jonkman (2002) Sidik, K. and Jonkman, J. N. (2002).  A simple confidence interval for meta-analysis.  Statistics in Medicine 21, 3153–3159.
  * Singh _and others_ (2005) Singh, K., Xie, M. and Strawderman, W.E. (2005).  Combining information from independent sources through confidence distributions.  The Annals of Statistics 33(1).
  * Sutton and Higgins (2008) Sutton, Alexander J. and Higgins, Julian P. T. (2008).  Recent developments in meta-analysis.  Statistics in Medicine 27(5), 625–650.
  * Veroniki _and others_ (2018) Veroniki, A. A., Jackson, D., Bender, R., Kuss, O., Langan, D., Higgins, J.P.T., Knapp, G. and Salanti, G. (2018).  Methods to calculate uncertainty in the estimated overall effect size from a random-effects meta-analysis.  Research Synthesis Methods 10(1), 23–43.
  * Xie _and others_ (2011) Xie, M., Singh, K. and Strawderman, W.E. (2011).  Confidence distributions and a unifying framework for meta-analysis.  Journal of the American Statistical Association 106(493), 320–333.


# A Functional approach for Two Way Dimension Reduction in Time Series

Aniruddha Rajendra Rao
Industrial AI Lab, Hitachi America, Ltd. R$\&$D
Santa Clara, CA
Haiyan Wang
Industrial AI Lab, Hitachi America, Ltd. R$\&$D
Santa Clara, CA
Chetan Gupta
Industrial AI Lab, Hitachi America, Ltd. R$\&$D
Santa Clara, CA

###### Abstract

The rise in data has led to the need for dimension reduction techniques,
especially in the area of non-scalar variables, including time series, natural
language processing, and computer vision. In this paper, we specifically
investigate dimension reduction for time series through functional data
analysis. Current methods for dimension reduction in functional data are
functional principal component analysis and functional autoencoders, which are
limited to linear mappings or scalar representations for the time series,
which is inefficient. In real data applications, the nature of the data is
much more complex. We propose a non-linear function-on-function approach,
which consists of a functional encoder and a functional decoder, that uses
continuous hidden layers consisting of continuous neurons to learn the
structure inherent in functional data, which addresses the aforementioned
concerns in the existing approaches. Our approach gives a low dimension latent
representation by reducing the number of functional features as well as the
timepoints at which the functions are observed. The effectiveness of the
proposed model is demonstrated through multiple simulations and real data
examples.

{Aniruddha.Rao, Haiyan.Wang<EMAIL_ADDRESS>

Keywords: Functional Data Analysis, Deep Learning, Autoencoder, Functional
Neural Network, Time series, Dimension Reduction.

##  1 Introduction

Nowadays, with the rapid advancements in information technology, data are
continuously collected at an astonishing frequency in many fields, such as
sensors installed on industrial equipment, economic factors, individual
health, and environmental exposures. The exponentially growing volume of data
poses challenges for storage, transfer, and analysis.

Dimension reduction plays a key role in solving this problem by reducing the
data in a meaningful manner. The objective of dimension reduction is to
mathematically reduce the dimensionality by projecting the data to a lower
dimensional subspace with minimum loss of information. The learned
mathematical mapping plays an effective role in not only sorting out which
variables are important but also how they interact with each other. Dimension
reduction also helps us to deal with the curse of dimensionality, ease the
transfer of information, reduce storage requirements, and reduce computation
time. The usefulness of dimension reduction makes it an active area of
research. Particularly, in the last few years, building dimension reduction
models have proven to be beneficial in many fields, like time series, text,
image, and video. This has attracted the interests of researchers from the
scientific community, especially where non-scalar types of data have become
prevalent [27, 13, 24, 8, 21]. Unfortunately, there is still a lot to be done
in dimension reduction for time series data.

The purpose of this paper is to consider the problem of dimension reduction in
time series data that occurs frequently in many areas of interest in industry
and research like economics, meteorology, finance, health, sensors, etc.
Examples include the weekly temperature, hourly sensor data, daily stock
returns, and monthly blood pressure readings of a patient. This is a difficult
and interesting problem as time series consists of random observations of the
same variables at different timestamps, and the information is not independent
nor linear because of the intricate temporal correlations. Mathematically, we
are building mapping from multiple chronologically measured numerical
variables within a certain time interval $\mathcal{S}$ to a few
chronologically measured numerical variables within the same time interval
$\mathcal{S}$.

In this paper, we consider a general setting similar to an autoencoder where
the input and output time series are the same and we learn the dimension
reduced form of the time series using the latent space representation. The
standard machine learning methods for dimension reduction, such as Principal
component analysis (PCA) and Autoencoders (AE), are designed to work with data
that are recorded at regularly spaced points and in the finite dimensional
space. These methods have the following limitations: ignoring intricate
temporal correlations, unable to capture complex relations, scalar
representation of time series data [2].

Functional data analysis (FDA) [15, 4, 9] considers this time series
information as functions over a continuum that are intrinsically infinite
dimensional. Functional principal component analysis (FPCA) [15, 4, 9] is the
functional counterpart to PCA from the multivariate setting. FPCA is a useful
dimension reduction tool that frequently serves as a key component in many
functional analysis [3, 5]. It allows us to do unsupervised learning of low-
dimensional representations of the time series data by considering the entire
time series as individual functional samples. FPCA helps to overcome some of
the limitations mentioned above but they learn a linear representation of the
data and therefore suffer from under fitting when the underlying mapping is
complex. There is also a sufficient dimension reduction branch in FDA [10, 11]
that learns a low-dimensional latent representation of the data which is non-
linear but they are supervised. This limits the applicability of these
methods. Also, they either don’t work or don’t perform well in the case of
multiple functional features.

Deep learning approaches have become very popular in the past few years. They
offer some of the state-of-the-art methods for learning nonlinear
representations of multiple types of data [23, 7]. More recently, deep
learning in functional data [18, 26, 25, 22, 17] has gained a lot of momentum.
Functional Autoencoders (FAE) [8] was inspired by Rossi et al. [19, 20],
introduced dimension reduction for functional data using deep learning. FAE
adapted AE to the functional setting. They showed that the FPCA is a special
case of FAE under certain conditions. While FAE has proven to beat a lot of
the state-of-the-art methods [8], unfortunately, it does not capture the true
nature of functional data as only some of the layers have weight functions and
the low dimensional latent representation of the time series data is scalar,
which is restrictive. So far, all the methods we discussed have some kind of
structural limitation, and to the best of our knowledge, the current works
tackle either decreasing the number of functional features or decreasing the
number of timepoints at which these curves are observed but there is very
limited work on solving both of them together.

In this paper, we formulate the dimension reduction for the time series
problem from the functional data analysis perspective. We innovatively
identify a general mathematical mapping between the functional features, based
on which a non-linear approach for two-way dimension reduction is proposed.
Specifically, we introduce Bi-Functional Autoencoders (BFAE), a novel
generalization of traditional autoencoders to functional data settings. The
contributions of this paper are summarized as follows:

  * •

We propose a general functional mapping from multivariate temporal variables
to themselves that embraces Functional Auto-encoders (FAE) as a special case.

  * •

We propose a new model that generalizes the well known neural network
autoencoders to the functional data setting while preserving the functional
nature of the data to address dimension reduction in time series.

  * •

We explain how our approach can capture non-linear relations while reducing
the dimensions in two-ways (number of features and time points) and,
therefore, called Bi-Functional Autoencoders (BFAE).

  * •

We demonstrate the effectiveness of the proposed approach in learning low
dimensional non-linear representation of the time series through simulation
experiments and two real-world examples.

##  2 Preliminaries

###  2.1 Notations and Prior Art

This paper aims to develop an approach to map and represent the multiple time
series variables (features) in a compact and efficient dimension reduced
manner by leveraging the temporal dependencies within and between the involved
variables. We begin by introducing useful notations, discussing the prior arts
and problem definition, before proceeding to give the details of Bi-Functional
Autoencoders (BFAE). Let us assume that for the $i^{th}$
($i\in\\{1,2,...,N\\}$) independent subject, we have $R$ features, that are
continuously recorded within a compact time interval
$\mathcal{S}\subseteq\mathbf{R}$. In particular, the observed timepoints of
the $r^{th}$ feature for the $i^{th}$ subject are given by a
$M^{(i,r)}_{s}$-dimensional vector
$\mathbf{S}^{(i,r)}=[S^{(i,r)}_{1},...,S^{(i,r)}_{j},...,S^{(i,r)}_{M^{(i,r)}_{s}}]^{T}$,
with $M^{(i,r)}_{s}$ representing the number of timepoints in the time series
and $S^{(i,r)}_{j}\in\mathcal{S}$ for
$i=1,...,n;r=1,...,R;j=1,...,M^{(i,r)}_{s}$. The corresponding time series are
denoted as
$\mathbf{X}^{(i,r)}=[X^{(i,r)}_{1},...,X^{(i,r)}_{j},...,X^{(i,r)}_{M^{(i,r)}_{s}}]^{T}$.
The subscript $M^{(i,r)}_{s}$ reflects the fact that the measuring timestamps
may vary across features and subjects. For the applicability of prior arts
like PCA and AE, the data needs to be regular time series, therefore, we
assume the number of timepoints across all features and subject is $M$. The
observed data can be denoted as
$\\{\mathbf{X}^{(i,1)},...,\mathbf{X}^{(i,R)}\\}_{i=1}^{N}$. We then
concatenate the $R$ temporal features as
$\mathbf{{X}}^{(i)}=[\mathbf{{X}}^{{(i,1)}^{T}},...,\mathbf{{X}}^{{(i,R)}^{T}}]^{T}$.
Given samples $\\{\mathbf{{X}}^{(i)}\\}_{i=1}^{N}$, the problem definition can
be stated as an unsupervised nonlinear multi-dimensional functional
representation learning problem, formally defined as follows:

$\mathbf{{Z}}^{(i)}=F(\mathbf{{X}}^{(i)})$ (1)

where
$\mathbf{Z}^{(i)}=\\{\mathbf{Z}^{(i,1)},...,\mathbf{Z}^{(i,R^{\prime})}\\}_{i=1}^{N}$
is a latent representation of the time series $\mathbf{{X}}^{(i)}$ with
$R^{\prime}<R$.

This problem formulation has a few disadvantages. Biases may get introduced
when performing data pre-processing for having the same number of timepoints
for each time series and these widely used prior arts have their own
limitations in solving the dimension reduction problems in time series data
when they assume them to be scalar values. These multivariate approaches are
unable to capture the temporal dependencies among $\mathbf{{X}}^{(i,r)}$,
$r=1,...,R$.

###  2.2 Functional Data Analysis and an alternate Formulation

In the paper so far, we have seen the dimension reduction problem from a
multivariate time series point of view. Let us look at an alternative problem
formulation using functional data analysis (FDA). In FDA, we learn from random
functions which are dynamically varying data over a continuum. There is a lot
of application for FDA in a wide variety of fields, like time series, sensor
data, image, and spatial data [14, 9]. We can use this ability of FDA for
dealing with continuous underlying curves $X^{(i,r)}(s),s\in\mathcal{S}$ that
is observed at discrete timepoints along a time interval. In functional
analysis, the input or output has to be functional curves if not both. For
given functional features
$\\{X^{(i,1)}(s),...,X^{(i,R)}(s),s\in\mathcal{S}\\}_{i=1}^{N}$, our problem
of finding a latent representation of the data can be given as follows:

$\mathbf{{Z}}^{(i)}(t)=F(\mathbf{{X}}^{(i)}(s)).$ (2)

where
$\mathbf{Z}^{(i)}(t)=\\{\mathbf{Z}^{(i,1)}(t),...,\mathbf{Z}^{(i,R^{\prime})}(t),t\in\mathcal{S}\\}_{i=1}^{N}$
is a latent representation of the time series $\mathbf{{X}}^{(i)(s)}$ with
$R^{\prime}<R$ and each time series is observed at $M^{\prime}$ timepoints
($M^{\prime}\leq M$).

The popularity of FDA in time series application is evident with the boom of
recent interest in this area of research. The assumption of smoothness can
deviate in FDA as shown in [28], allowing us to directly apply different FDA
techniques to time series data. Also, compared to the common sequential models
(RNNs, LSTMs) which learn fixed parameters over time, FDA enables us to learn
feature effects that change over the interval $\mathcal{S}$ and build more
efficient models [25]. In the functional field, we have FPCA [6, 1] which can
encode the information in a low-dimensional latent space but the mapping is
linear. This can be a limiting factor because, in many real-world
applications, the multiple functional features can have complex relations not
just among them but across different timepoints as well.

The growing interest in deep learning and FDA has resulted in multiple
successive deep functional models [19, 26, 16] for the different analytical
tasks. FAE [8] is a deep learning approach for dimension reduction, a
generalization of the autoencoder from the vector space to the functional
setting. In this network, only the first and last layer can accommodate
functions with the help of functional neurons given by [19, 20], and the rest
of the network consists of traditional neurons. Because of this, FAE can only
give a latent representation of functional data in the scalar form, which is
limiting. It is very common for functional features to have complex relations
that can’t be embedded in a scalar form. Thus, the scalar representation of
functional data is inadequate.

In the next section, we propose a non-linear function-on-function autoencoder
that leverages the power of fully connected Neural Networks while keeping the
functional identity of the time series.

##  3 Proposed Bi-Functional Autoencoder Model

###  3.1 Our approach

Figure 1: General architecture for our proposed Bi-Functional Autoencoder

The main goal of this paper is to learn a mapping $\mathbf{f}(\cdot)$ for the
time series data to a low dimensional latent space as shown in Figure 1. To
achieve this goal we make use of the idea of an autoencoder and Functional
Neural Networks [19, 16, 17]. An autoencoder (AE) is a popular dimension
reduction technique which can learn non-linear latent presentation in a
$R^{\prime}$-dimensional space $\mathbf{R}^{R^{\prime}}$ from a
$R$-dimensional vector-valued input space $\mathbf{R}^{R}$ (where
$R^{\prime}<R$) using the help of an encoder. The decoder part of the AE helps
to map the learned latent representation back to the $R$-dimensional vector-
valued input space $\mathbf{R}^{R}$. The input and output of an AE are the
same and we learn by measuring the reconstruction error of the output compared
to the input. We want to use this idea but in a functional setting where we
have an encoder that is mapping from a $R$-dimensional functional space to a
$R^{\prime}$-dimensional functional space and the decoder maps
$R^{\prime}$-dimensional functional space back to the $R$-dimensional
functional space. We make use of continuous neurons defined in [16, 17] to
enable us to achieve this task in a functional setting. The continuous neurons
take in functional inputs and produce a functional output, thus leading to the
preservation of the functional nature of the data.

The framework for our Bi-Functional Autoencoder (BFAE) Model, as seen in
Figure 1, involves development of a continuous mapping from layer to layer by
using multiple continuous hidden layers consisting of multiple continuous
neurons as defined in [16, 17]. In Figure 1, we have three types of layers: an
input layer, continuous hidden layers, and a continuous output layer. The
input layer is the first layer that takes in the functional features, the
continuous hidden layers consist of multiple continuous neurons, and the
continuous output layer provides the reconstructed input functional features
using the continuous neurons. The model learns with the help of functional
weights that are continuous over time (or some other continuum). These
functional weights consists of both univariate and bivariate functions, that
is, we learn functions as well as surfaces. The key idea here is that we want
to preserve the functional nature of the data throughout the whole
reconstruction process across the network. This will give a richer structure
to transform the data into low dimensional latent representation while
exploiting the domain information and continuity of the functional features.
We define the $l^{th}$ continuous hidden layer and its $r^{th}$ continuous
neuron as:

$\begin{split}H^{(i,r)}_{(l)}(s)&=\sigma\Big{(}b^{(r)}_{(l)}(s)+\sum_{j=1}^{J}\int
w^{(r,j)}_{(l)}(s,t)H^{(i,j)}_{(l-1)}(t)dt\Big{)}\end{split}$ (3)

where $l=1,2,3,...,L$, $H^{(i,r)}_{(0)}(s)=X^{(i,r)}(s)$,
$H^{(i,r)}_{(L)}(s)=\widehat{X^{(i,r)}(s)}$,
$b^{(r)}_{(l)}(s)\in\mathcal{L}^{2}(\mathcal{S})$ is the unknown univariate
intercept function,
$w^{(r,j)}_{(l)}\in\mathcal{L}^{2}(\mathcal{S}\times\mathcal{S})$ is the
bivariate parameter function for the $r^{th}$ continuous neuron in the
$l^{th}$ continuous hidden layer coming from the $j^{th}$ continuous neuron of
the $(l-1)^{th}$ continuous hidden layer and $\sigma(\cdot)$ is a standard
activation function. The value of $J$, that can be considered as the number of
features in each continuous hidden layer, is fixed for the input layer and the
continuous output layer as $R$. We can tune $J$ to any scalar value for the
rest of the continuous hidden layers. Similarly, the timepoints at which the
output of the continuous neurons is observed is fixed for the input layer and
the continuous output layer and can be changed for all the other continuous
neurons in the remaining continuous hidden layers.

We will now look into the working of our network. We use the continuous
neurons defined by Equation 3, to get from the input layer to a dimension
reduced form of the input at an intermediate continuous hidden layer (say, it
is the $l^{\prime th}$ layer). This is the encoder part of the BFAE approach
called as the functional encoder. We can get the low dimensional
representation of the inputs in the $l^{\prime th}$ continuous hidden layers
as $\mathbf{Z}^{(i)}(s)=H^{(i)}_{(l^{\prime})}(s)$ where
$\mathbf{Z}^{(i)}(s)=\\{\mathbf{Z}^{(i,1)}(s),...,\mathbf{Z}^{(i,R^{\prime})}(s),s\in\mathcal{S}\\}_{i=1}^{N}$
and
$H^{(i)}_{(l^{\prime})}(s)=\\{H^{(i,1)}_{(l^{\prime})}(s),...,H^{(i,R^{\prime})}_{(l^{\prime})}(s),s\in\mathcal{S}\\}_{i=1}^{N}$,
$R^{\prime}<R$ and the function $\mathbf{Z}^{(i)}$ is observed at
$M^{\prime}\leq M$ timepoints. We move from this latent representation of the
input towards the continuous output layer using the continuous neurons in the
continuous hidden layers to get the reconstructed functional input values.
This part of our approach is called the functional decoder. Therefore, we have
developed a mapping from a $R$-dimensional functional space to
$R^{\prime}$-dimensional functional space using the functional encoder, and
the functional decoder maps $R^{\prime}$-dimensional functional space back to
the $R$-dimensional functional space. This latent representation
$\mathbf{Z}^{(i)}(s)$ is very flexible in nature, where we can set the number
of functional features as $R^{\prime}$ and the number of timepoints observed
as $M^{\prime}$, according to the task in hand.

For simplicity, let us assume that Figure 1 has $R$ functional features in the
input layer, $L$ continuous hidden layers, each with $J$ incoming and $R$
outgoing connections. BFAE is similar to an autoencoder in nature, but for
functional data, we have to define the number of continuous hidden layers, the
number of continuous neurons in each of the continuous hidden layers, and the
number of timepoints at which these functions are observed for each continuous
neuron. The activation function used in the continuous neurons is ReLU, tanh,
sigmoid/logistic, or linear (continuous output layer). The forward propagation
of the network is straightforward with the help of Equation 3. In the
backpropagation phase, we learn using the functional surfaces (bivariate) and
the functional intercepts (univariate) directly. We need the functional
gradients to learn the functional parameters of the BFAE network. These
functional gradients measure the change in a functional to a change in a
function on which the functional depends.

We use the functional gradients given below to optimize our network’s
functional parameters. The optimization approach we use is the traditional
gradient descent. The necessary assumptions and mathematical tools are
discussed in [19, 26, 12]. We use Fréchet derivatives, from the calculus of
variation, for computing the functional gradients. The partial derivatives
needed for the backpropagation are as follows:

$\begin{split}\frac{\partial H^{(i,r)}_{(l)}}{\partial
b^{(r)}_{(l)}}(s)&=\sigma^{\prime}\Big{(}b^{(r)}_{(k)}(s)+\sum_{j=1}^{J}\int
w^{(r,j)}_{(l)}(s,t)H^{(i,j)}_{(l-1)}(t)dt\Big{)}\\\ \frac{\partial
H^{(i,r)}_{(l)}}{\partial
w^{(r,j)}_{(l)}}(s,t)&=\sigma^{\prime}\Big{(}b^{(r)}_{(l)}(s)+\sum_{j^{\prime}=1}^{J}\int
w^{(r,j^{\prime})}_{(l)}(s,t^{\prime})H^{(i,j^{\prime})}_{(l-1)}(t^{\prime})dt^{\prime}\Big{)}\\\
&\times\int\frac{\partial}{\partial
w^{(r,j)}_{(l)}(s,t)}w^{(r,j)}_{(l)}(s,t)H^{(i,j)}_{(l-1)}(t)dt\\\
&=\sigma^{\prime}\Big{(}b^{(r)}_{(l)}(s)+\sum_{j^{\prime}=1}^{J}\int
w^{(r,j^{\prime})}_{(l)}(s,t^{\prime})H^{(i,j^{\prime})}_{(l-1)}(t^{\prime})dt^{\prime}\Big{)}\\\
&\times H^{(i,j)}_{(l-1)}(t)\\\ \end{split}$ (4)

where $l=1,2,3,...,L$, $H^{(i,r)}_{(0)}(s)=X^{(i,r)}(s)$,
$H^{(i,r)}_{(L)}(s)=\widehat{X^{(i,r)}(s)}$ and $\sigma^{\prime}(\cdot)$
represents the first derivative of $\sigma(\cdot)$. In the backpropagation
phase, we pass through the network backward from the continuous output layer
to the input layer and calculate the partial derivatives of the continuous
neurons with respect to the bivariate functional weight and the functional
intercept. The loss function that we minimize is $\mathcal{L}(\theta)$, where
$\theta$ is the collection of all functional parameters in the BFAE network.
The loss function is given as follows:

$\mathcal{L}(\theta)={\frac{1}{N}\sum_{i=1}^{N}\sum_{r=1}^{R}\int\Big{(}X^{(i,r)}(s)-\widehat{X^{(i,r)}(s)}\Big{)}^{2}dt}$
(5)

The number of continuous hidden layers, the number of continuous neurons in
each of the continuous hidden layers, and the number of observed timepoints
for each continuous neuron can be considered as hyperparameters. We can also
adjust our approach to accommodate irregular functional data or scalar
variables. All these offer a lot of flexibility for our approach to work with
different kinds of problems and adjust things according to the downstream
process for different analytical tasks like prediction, classification,
clustering, forecasting, and more.

###  3.2 Connection of our approach to prior arts

Some of the current dimension reduction methods can be represented as a
special case of BFAE. As we know, a single layer autoencoder (AE) with a
linear activation function is analogous to PCA. In the same way, we can adjust
the parameters of our model to act like FPCA. The FPCA is a special case of
BFAE when our model has the following specifications: linear activation
functions, a single continuous hidden layer, and the functional weights are
constrained to be orthonormal. The Functional Autoencoder (FAE) [8] paper
discusses how AE is a special case of FAE as FAE replaces the scalar weights
and inner products of the AE with functional weights and inner products. In
the same manner, FAE is a special case of our approach, BFAE. If we specify
the functional encoder in our model to learn a scalar ($M^{\prime}=1$) latent
representation of the functional features, our approach essentially acts like
FAE. Therefore, our approach is a generalization of most of the existing
methods like PCA, FPCA, AE, and FAE.

##  4 Numerical Experiments

In this section, we proceed to apply the proposed model to multiple simulation
scenarios and two real-world problems. We compare our results against several
state-of-the-art approaches, like PCA, AE, and FPCA, to show the effectiveness
of BFAE. We demonstrate the following objectives through our results: 1) Our
approach is effectively trained by the derived functional gradients 2) BFAE
enables learning of the relation between the functional features 3) We show
the efficacy of our approach by outperforming other methods.

###  4.1 Simulations

In the first experiment, we consider an individual predictor function (i.e.
$R=1$) which is observed on a dense and regular grid. We generate $N=100,1000$
iid random curves $\left\\{X^{(1)}(t),\cdots,X^{(n)}(t)\right\\}$ from a
Gaussian process with mean $0$ and covariance given as follows:

$C_{X}(t,s)=\frac{\sigma^{2}}{\Gamma(\nu)2^{\nu-1}}\left(\frac{\sqrt{2\nu}|t-s|}{\rho}\right)^{\nu}K_{\nu}\left(\frac{\sqrt{2\nu}|t-s|}{\rho}\right).$
(6)

where this is the Matérn covariance function and $K_{\nu}$ is the modified
Bessel function of the second kind. The value of $\rho=0.5,$ $\nu=5/2$, and
$\sigma^{2}=1.$ These curves are realized at equally-spaced timepoint on the
interval $[0,1]$ where $M=50,250$ timepoints. We have specifically chosen a
case where $N<M$ to check if it causes any issues with any of the approaches.

After generating the curves using a Gaussian process, we add some random noise
($\epsilon_{i}(t)$) to them. The results are averaged over 100 simulations,
and we divide the 100 samples into 80 for training and 20 for testing. The
same ratio split is used for the case of 1000 samples. We report the Root Mean
Square Error (RMSE) for the reconstruction of the functional features given by
Equation (7) and in the case of the scalar approaches, we just report the
standard RMSE ignoring the temporal aspect. For our approach in the first
experiment, we consider multiple combinations of continuous hidden layers
($L$= 1, 3) and only one continuous neuron ($J=1$). The grid (s) value for
latent representation in the continuous neurons is either $M$ or
$M^{\prime}=M/5$. We keep the network of AE similar to BFAE, and for PCA and
FPCA we follow the $99\%$ variance of explained rule. Note that we can select
$M^{\prime}$ to be any scalar value in the range of ($1,M$) but we select this
particular value ($M^{\prime}=M/5$) to show the flexibility of our approach.

$RMSE=\sqrt{\frac{1}{N}\sum_{i=1}^{N}\sum_{r=1}^{R}\int\Big{(}X^{(i,r)}(s)-\widehat{X^{(i,r)}(s)}\Big{)}^{2}dt}$
(7)

M | 50 | 250 | 50 | 250
---|---|---|---|---
Method/N | 100 | 1000
PCA | 0.483 | 0.418 | 0.469 | 0.398
AE | 0.409 | 0.443 | 0.390 | 0.406
FPCA | 0.390 | 0.351 | 0.356 | 0.343
BFAE | 0.301 | 0.265 | 0.102 | 0.098
BFAE (M’) | 0.331 | 0.274 | 0.117 | 0.104
Table 1: Comparing RMSE of different methods for capturing a single time
series function.

We observe from Table I that PCA and AE perform similarly on the temporal
curves. There is some gain in performance for FPCA but our approach performs
the best irrespective of the sample size ($N$) and the number of timepoints
($M$). The difference in performance increases when the sample size and the
number of timepoints increases, as our approach is a deep learning model, it
learns better with more information. When we set the number of observed
timepoints as $M^{\prime}$, our approach loses some performance because we are
restricting the true data to a smaller latent space. But, even under this
scenario of representing the true curves using 10 or 50 timepoints rather than
50 or 250, we are still performing better than all the other approaches.

Figure 2: Curve reconstruction comparison for different methods when $N=1000$,
$R=1$, and $M=50$

Figure 2 shows us the reconstructed curves for a random sample from different
approaches for the simulation setting of $N=1000$, $R=1$, and $M=50$. We see a
similar story compared to Table I where PCA is performing poorly, and the
shape is very different from the truth. PCA is a linear approach and can’t
capture the temporal information in the data. AE and FPCA perform better than
PCA by getting the general shape correct but are still far off from the truth.
Our approach, BFAE, irrespective of the number of timepoints observed, is able
to reconstruct the truth well and captures the shape of the curve correctly
over the whole time interval.

We extend the first experiment to include more functional features. In the
second experiment, we increase the number of features to $R=10$. We follow the
same procedure to generate these 10 curves for the same set of values of $N$
and $M$. The results discussed below are an average of 100 simulations after
adding random noise and using the same splits for training and testing. For
our approach, we consider multiple combinations of continuous hidden layers
($L$= $1$, 3), continuous neurons ($J$= 1, 2, 4) and we set $R^{\prime}=4$.
The grid (s) value for latent representation in the continuous neurons is
either $M$ or $M^{\prime}=M/5$.

M | 50 | 250 | 50 | 250
---|---|---|---|---
Method/N | 100 | 1000
PCA | 1.527 | 1.525 | 1.518 | 1.518
AE | 0.471 | 0.490 | 0.489 | 0.484
FPCA | 0.386 | 0.417 | 0.365 | 0.404
BFAE | 0.328 | 0.292 | 0.109 | 0.132
BFAE (M’) | 0.359 | 0.314 | 0.147 | 0.169
Table 2: Comparing RMSE of different methods for capturing multiple ($R=10$)
time series functions.

Table II shows that the increase in the functional features results in a
deterioration in performance for PCA. While the other two approaches perform
better than PCA, our approach performs the best irrespective of the sample
size ($N$) and the number of timepoints ($M$). The difference in performance
again increases when the sample size and number of timepoints increases. We
observe that the errors have increased marginally compared to the previous
table as we have added more features. The behavior of decreasing the number of
timepoints to $M^{\prime}$ is similar as well.

###  4.2 Real Data applications

We consider two real data sets. The first is the speech recognition data from
TIMIT (available at http://statweb.stanford.edu/tibs/Ele-mStatLearn/) where we
have speech signals for different phonemes. Audio information is available
abundantly, but they are measured at a very high frequency requiring high
storage capacity. A low dimensional latent representation of this data would
be very useful in the industry. The other example deals with the relation
between electricity demand and temperature in the city of Adelaide, Australia.
This data is important because of the high costs related to the storage of
electricity and understanding the effect of temperature on demand can lead to
information based operational steps.

For speech recognition data, we set up the experiment similar to [10, 16],
where we have voice signals ($R$=1) for two phonemes (response) transcribed as
follows: “aa” as the vowel in “dark” and “ao” as the first vowel in “water”.
We build a classifier model with the help of FLM after applying dimension
reduction to the voice signals. Figure 3 shows that it is difficult to
separate the two groups where we see the two phoneme curves computed by a log-
periodogram of 150 ($M=150$). Each phonemes information is recorded over 150
points, reduction of such data is beneficial, especially with the potential to
expand this low dimension representation to much higher voice signals. We have
800 functional samples, we split these samples into 640 samples for training
and 160 samples for testing.

The Adelaide data (available in $fds$ package in R-software) has the
temperature and electricity demand records from 7/6/1997 to 3/31/2007 (508
weeks) measured daily at half-hourly rate ($M=48$). We consider each day of
the week as a feature ($R$=7) and try to map the relation between temperature
and demand using Functional Linear model (FLM) [4, 9, 14]. Before we use FDA
to model this relation, we use BFAE to represent the temperature curves in a
compact low dimensional form. We split this data into 400 samples for training
and 108 samples for testing. Figure 4 shows the half-hourly temperature and
electricity demand (Megawatts) for the whole 508 weeks. We can see from the
image that the temperature curve for each day of the week follows a similar
pattern and learning that pattern can lead to a reduction in dimension with
minimum loss of information.

Method/Data | Phonemes | Adelaide
---|---|---
PCA | 14.960 | 18.633
AE | 1.790 | 2.354
FPCA | 1.242 | 0.636
BFAE | 1.126 | 0.581
BFAE (M’) | 1.128 | 0.615
Table 3: Comparing RMSE of different methods for capturing prediction
functions for the real data.

We can observe from Table III that PCA again is having difficulties in
reducing the information for both the real data sets. While AE and FPCA are
performing better than PCA, BFAE performs the best. For phonemes data, we even
reduced the points to $M^{\prime}=30$ and still BFAE performs better than
other approaches. The difference between BFAE with timepoints as $M$ and
$M^{\prime}$ is negligible, indicating that the low dimension representation
is very rich in information. While for Adelaide data, we reduce the number of
features to 4 ($R^{\prime}=4$) and the timepoints to $M^{\prime}=12$ and
produce the best results using our approach with timepoints as $M$. Our
intuition to reduce the temperature feature is certainly validated as seen
from Figure 4 and Table III.

Modeling results can be seen in Table IV where we build a classifier for the
phonemes data, to predict if the voice signal is for “aa” or “ao”. For the
Adelaide data, we map the temperature information for the 7 days of the week
to predict the electricity demand for the 7 days of that same week. Both these
models are built with the help of the original data and the low dimension
latent representation learned using BFAE. We ignore the other approaches as
the goal here is to demonstrate that our approach gives modeling results at
least as good as the original data by capturing the signals. We can see from
Table IV that our approach not only does a good job of modeling but also
performs as well as the original data and has a lower tendency to over fit as
it contains less noise. The error increases in both cases as we reduce the
timepoints to $M^{\prime}$ but the performance is still competitive.

Real Data | Data | Original | BFAE | BFAE ($M^{\prime}$)
---|---|---|---|---
| train | 0.175 | 0.185 | 0.205
Phonemes | test | 0.200 | 0.190 | 0.210
| train | 188.363 | 164.105 | 199.758
Adelaide | test | 220.739 | 184.518 | 215.294
Table 4: Comparison of errors (Classification error for Phonemes and RMSE for
Adelaide) of different methods for modeling the response.

Overall, our approach does the best job in reducing the information into a low
dimensional latent space while still maintaining a competitive performance at
modeling different tasks compared to the original data. We are able to
represent the information in a compact data rich manner and are also able to
reconstruct it back to the original scale without much loss of information.

Figure 3:  Voice signal curves for the Phonemes data

Figure 4:  Day curves (Monday to Sunday, top to Bottom) of temperature and
electricity demand for Adelaide

##  5 Conclusions

In this paper, we proposed a novel model for multivariate time series
dimension reduction. Our approach helps to perform two-way dimension reduction
by reducing the number of features and the number of timepoints at which the
time series is observed. The proposed Bi-Functional Autoencoder (BFAE) reduces
the input into a low dimension latent representation using a functional
encoder and reconstructs the information back using a functional decoder. Our
proposed approach has a lot of advantages compared to current methods,
including its ability to represent the time series data in a flexible manner,
capture timely varying correlations among features, ability to deal with
different kinds of data, and capture non linear relations. Along with superior
simulation results, we showed the proficiency of our approach in two real-
world examples in comparison with several common practices in the prior art.
Our approach produced smaller errors in reconstruction and modeled the data as
well as the original information. We expect the proposed model to be widely
used in diverse real-world problems where the goals are to reduce the transfer
load of huge amount of time series data, store a large amount of time series
information effectively, and reduce the computation cost of different
analytical tasks while retaining the performance.

## References

  * [1] Daniel Backenroth, Jeff Goldsmith, Michelle D. Harran, Juan C. Cortes, John W. Krakauer, and Tomoko Kitago.  Modeling motor learning using heteroscedastic functional principal components analysis.  Journal of the American Statistical Association, 113(523):1003–1015, 2018.  PMID: 30416231.
  * [2] Wei Chen, Hythem Sidky, and Andrew L. Ferguson.  Capabilities and limitations of time-lagged autoencoders for slow mode discovery in dynamical systems.  The Journal of Chemical Physics, 151(6):064123, 2019.
  * [3] Jeng-Min Chiou, Yu-Ting Chen, and Ya-Fang Yang.  Multivariate functional principal component analysis: A normalization approach.  Statistica Sinica, pages 1571–1596, 2014.
  * [4] F. Ferraty and Y. Romain.  The Oxford Handbook of Functional Data Analysis.  Oxford University Press, 01 2011.
  * [5] Clara Happ and Sonja Greven.  Multivariate functional principal component analysis for data observed on different (dimensional) domains.  Journal of the American Statistical Association, 113(522):649–659, 2018.
  * [6] Clara Happ and Sonja Greven.  Multivariate functional principal component analysis for data observed on different (dimensional) domains.  Journal of the American Statistical Association, 113(522):649–659, 2018.
  * [7] G. E. Hinton and R. R. Salakhutdinov.  Reducing the dimensionality of data with neural networks.  Science, 313(5786):504–507, 2006.
  * [8] Tsung-Yu Hsieh, Yiwei Sun, Suhang Wang, and Vasant Honavar.  Functional Autoencoders for Functional Data Representation Learning, pages 666–674.  arxiv, 2021.
  * [9] Piotr Kokoszka and Matthew Reimherr.  Introduction to Functional Data Analysis.  New York: Chapman and Hall/CRC, 2018.
  * [10] Bing Li and Jun Song.  Nonlinear sufficient dimension reduction for functional data.  Annals of Statistics, 45(3):1059–1095, June 2017.  Copyright: Copyright 2017 Elsevier B.V., All rights reserved.
  * [11] Bing Li and Jun Song.  Dimension reduction for functional data based on weak conditional moments.  The Annals of Statistics, 50(1):107 – 128, 2022.
  * [12] Peter J. Olver.  Introduction to the calculus of variations.  In Introduction to the Calculus of Variations, 2014.
  * [13] Stefan Petscharnig, Mathias Lux, and Savvas Chatzichristofis.  Dimensionality reduction for image features using deep learning and autoencoders.  In Proceedings of the 15th International Workshop on Content-Based Multimedia Indexing, CBMI ’17, New York, NY, USA, 2017. Association for Computing Machinery.
  * [14] James O Ramsay.  Functional data analysis.  Wiley Online Library, 2006.
  * [15] J.O. Ramsay and B.W. Silverman.  Functional Data Analysis.  Springer series in statistics. Springer, 1997.
  * [16] Aniruddha Rajendra Rao and Matthew Reimherr.  Modern non-linear function-on-function regression, 2021.
  * [17] Aniruddha Rajendra Rao and Matthew Reimherr.  Non-linear functional modeling using neural networks, 2021.
  * [18] Aniruddha Rajendra Rao, Qiyao Wang, Haiyan Wang, Hamed Khorasgani, and Chetan Gupta.  Spatio-temporal functional neural networks.  arXiv preprint arXiv:2009.05665, 2020.
  * [19] Fabrice Rossi, Brieuc Conan-Guez, and François Fleuret.  Functional data analysis with multi layer perceptrons.  In Proceedings of IJCNN, pages 2843–2848. Citeseer, 2002.
  * [20] Fabrice Rossi, Nicolas Delannay, Brieuc Conan-Guez, and Michel Verleysen.  Representation of functional data in neural networks.  Neurocomputing, 64:183–210, 2005.
  * [21] Nitish Srivastava, Elman Mansimov, and Ruslan Salakhutdinov.  Unsupervised learning of video representations using lstms.  In Proceedings of the 32nd International Conference on International Conference on Machine Learning - Volume 37, ICML’15, page 843–852. JMLR.org, 2015.
  * [22] Barinder Thind, Kevin Multani, and Jiguo Cao.  Deep learning with functional inputs, 2020.
  * [23] Michael Tschannen, Olivier Frederic Bachem, and Mario Lučić.  Recent advances in autoencoder-based representation learning.  In Bayesian Deep Learning Workshop, NeurIPS, 2018.
  * [24] Laurens van der Maaten, Eric Postma, and H. Herik.  Dimensionality reduction: A comparative review.  Journal of Machine Learning Research - JMLR, 10, 01 2007.
  * [25] Qiyao Wang, Haiyan Wang, Chetan Gupta, Aniruddha Rajendra Rao, and Hamed Khorasgani.  A non-linear function-on-function model for regression with time series data.  In 2020 IEEE International Conference on Big Data (Big Data), pages 232–239, 2020.
  * [26] Qiyao Wang, Shuai Zheng, Ahmed Farahat, Susumu Serita, Takashi Saeki, and Chetan Gupta.  Multilayer perceptron for sparse functional data.  In 2019 International Joint Conference on Neural Networks (IJCNN), pages 1–10. IEEE, 2019.
  * [27] Yasi Wang, Hongxun Yao, and Sicheng Zhao.  Auto-encoder based dimensionality reduction.  Neurocomputing, 184:232–242, 2016.  RoLoD: Robust Local Descriptors for Computer Vision 2014.
  * [28] 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.