∎

11institutetext: Yiyan Yao 22institutetext: School of Information and
Communication Engineering, Beijing University of Posts and Telecommunications,
P. O. Box 101, Xitucheng Road No. 10, Haidian District, Beijing China 100876
22email<EMAIL_ADDRESS>33institutetext: Xin-long Luo 44institutetext:
Corresponding author. School of Information and Communication Engineering,
Beijing University of Posts and Telecommunications, P. O. Box 101, Xitucheng
Road No. 10, Haidian District, Beijing China 100876
Tel.: +086-13641229062
44email<EMAIL_ADDRESS>
# Improving Vertical Positioning Accuracy with the Weighted Multinomial
Logistic Regression Classifier

Yiyan Yao  Xin-long Luo $\ast$

(Received: date / Accepted: date)

###### Abstract

In this paper, a method of improving vertical positioning accuracy with the
Global Positioning System (GPS) information and barometric pressure values is
proposed. Firstly, we clear null values for the raw data collected in various
environments, and use the 3$\sigma$-rule to identify outliers. Secondly, the
Weighted Multinomial Logistic Regression (WMLR) classifier is trained to
obtain the predicted altitude of outliers. Finally, in order to verify its
effect, we compare the Multinomial Logistic Regression (MLR) method, the WMLR
method, and the Support Vector Machine (SVM) method for the cleaned dataset
which is regarded as the test baseline. The numerical results show that the
vertical positioning accuracy is improved from 5.9 meters (the MLR method),
5.4 meters (the SVM method) to 5 meters (the WMLR method) for 67% test points.

###### Keywords:

vertical positioning data correction parameter estimation multinomial logistic
regression support vector machine global positioning system

##  1 Introduction

In recent years, the performance of the Global Positioning System (GPS) is
excellent in outdoor environments RGG2015 . When users are outdoors, their
locations can be obtained accurately through GPS. However, the GPS signals are
blocked by the buildings and other obstacles, which result in large indoor
positioning errors. Thus, the indoor positioning accuracy is often challenged,
especially in the vertical direction. In the meantime, the space that we are
living in is filled with many high-rise buildings and our most activities are
indoors. Considering the practical requirement and the poor indoor positioning
performance, researchers have tried many methods to improve the vertical
positioning accuracy, such as the WiFi-based localization technology
DZYXKY2018 ; LLH2017 ; ZHLJX2016  and the barometer-based positioning
technology XWJC2015 .

On the other hand, the GPS chip has been embedded in the most mobile
terminals, which provides the location and timing information such as time,
latitude, longitude, speed and altitude. Therefore, based on the GPS
information, many researchers put forward some effective methods to improve
the positioning accuracy of the low-cost GPS about 4 meters to 10 meters in
several experiments IK2014 . Huang and Tsai propose an approach to calibrate
the GPS position by using the context awareness technique from the pervasive
computing and improve the positioning accuracy of GPS effectively HT2008 . The
machine learning techniques are applied to assess and improve the GPS
positioning accuracy under the forest canopy in ORMMS2011 .

In this paper, we provide another machine learning technique ALTMY2018 ;
AASMGK2019 ; AZISYDML2019  based on the Multinomial Logistic Regression (MLR)
method KS2016 ; MGB2008  for the vertical positioning problem. The research
data are measured by many different user equipments and provided by Huawei
Technologies Company, some data of which include the GPS three-dimensional
information and the barometric pressure values, and Some data of which miss
the GPS information or the barometric pressure values. We preprocess the
research data firstly. Consequently, we identify the abnormal data with the
$3\sigma$-rule and clear them. Meanwhile, some noises arise from the
inaccurate data records and the different reference standards of different
kinds of user equipments. These intrinsic noises lead to the poor distribution
law between the air pressure and the corresponding altitude. In order to
overcome these noise effects, we convert this vertical positioning problem
into a classification problem and revise the weighted MLR method to improve
its vertical positioning accuracy. Finally, in order to verify the effect of
the Weighted Multinomial Logistic Regression (WMLR) method, we compare the MLR
method, the WMLR method, and the Support Vector Machine (SVM) method CL2001 ;
CL2013 ; CTS2017  for this vertical positioning problem. The numerical results
show that the vertical positioning accuracy of the cleaned data is improved
from 5.9 meters (the MLR method), 5.4 meters (the SVM method) to 5 meters (the
WMLR method) for 67% test points.

The rest of the paper is organized as follows. In section 2, some related
works are discussed. In section 3, we describes the methodology of the data
cleaning, the outlier detection and the data correction based on the WMLR
classifier. In section 4, we describe the data source and compare the MLR
method, the WMLR method and the SVM method for the cleaned data which is
regarded as the test baseline. The promising numerical results are also
reported. Finally, some conclusions and the further works are discussed in
section 5.

##  2 Related works

In the field of improving the indoor vertical positioning accuracy, many
studies have been conducted. The related works can be roughly divided into two
categories: the Received Signal Strength Strength Indication (RSSI) based
methods and the barometric pressure based methods.

The RSSI of the Wi-Fi and the cellular network based methods use the collected
the RSSI and build the database of the fingerprints for the floor positioning
BSHL2016 ; WLSL2014 ; ZLC2012 . Some researchers consider the locations of the
Wi-Fi access points to determine the floor GBRB2014 . In BSHL2016 , the
experimental data are collected from one or two buildings and the collecting
device is fixed. They use the collected RSSI information and the pressure data
to estimate the floor. In those papers, since the RSSI information is local,
when the experimental environment changes, the training data need to be
collected by hand and the discriminant parameters need to be trained again.

Since there are many Wi-Fi access points distributed in a crowded indoor
environment and the wall cannot completely obstruct the signals, the signal
interference and fluctuation of different floors will result in the inaccurate
estimation. Some researchers propose the barometric altimetry for the floor
determination. In XWJC2015 , Xia et al. give a method based on the multiple
reference barometers for the floor positioning in buildings and their method
can give an accurate floor level. The disadvantages of their method are that
the height thresholds should be given in the floor determination and they are
sensitive to the local pressure conditions.

In CTS2017 , Chriki et al. use the SVM method based on the RSSI measurements
for the zoning localization problem. In ALTMY2018 , Adege et al. propose an
outdoor and indoor positioning method based on the hybrid of SVM and deep
neural network algorithms according to the RSSI of the Wi-Fi. Since the SVM
method only considers the support vector and the few points which are most
relevant are used to make the classification, its classification result may be
ineffective when the level of noise is high. The positioning method based on
the deep neural network ALTMY2018 ; HLL2017  requires a very large amount of
data to perform better than other techniques, and it requires expensive GPUs
and multiple devices to train complex models. The MLR method considers all
training data points which smooth the noise such that the MLR method can
handle the high level of noise of the training data. Furthermore, the MLR
method can be used to handle the large scale problem K2019 . Therefore, in
consideration of the performance gain of the weighted positioning algorithm
LLH2017 , we choose the MLR method with the weighted technique as the vertical
positioning method based on the GPS and barometric pressure information of the
user equipments.

##  3 The methodology

Our positioning method is composed of several stages, including the data
cleaning, the outlier detection, the data correction and the prediction of
vertical altitude for the test feature vector. We described these procedures
in the following subsections.

###  3.1 Data cleaning

The raw dataset is measured at different places with different user
equipments. In the dataset, many data miss the air pressure values due to some
mobile devices without the barometers. We delete these data of the missing air
pressure values firstly. Additionally, there are some abnormal data which
deviate too far from the average value of the dataset and it is shown as
follows. Assume that an average sea level pressure is 1013.25 hPa and the
corresponding temperature is 15∘C, then the air pressure value and its
corresponding altitude have the following relationship ZF2014 :

$\displaystyle h=44330.8-4946.54p^{0.1902632},$ (1)

where the unit of altitude $h$ is meter, and the unit of the air pressure
value $p$ is hpa. From formula (1), it is not difficult to find that the
barometric pressure value and the corresponding altitude are the inverse
relationship. However, from Fig. 1, we find that the distribution between the
air pressure values and the corresponding altitudes of the given data is
irregular. Therefore, we conclude that there exists the data drift in the
given real test data. Thus, we use the 3$\sigma$-rule to exclude the abnormal
data as follows W2004 :

$\displaystyle X\;\text{is thrown away when}\;|X-\mu|\geq 3\sigma,$

where the mean $\mu$ and the standard deviation $\sigma$ are computed by the
following formula:

$\displaystyle\mu=\frac{1}{n}\sum_{i=1}^{n}X_{i}\;\text{and}\;\sigma=\sqrt{\frac{1}{n-1}\sum_{i=1}^{n}(X_{i}-\mu)^{2}}.$

After performing the 3$\sigma$-rule, we eliminate the large deviation data and
the $99.73\%$ data are retained.

Figure 1: The distribution of the pressure values and the corresponding
altitudes.

###  3.2 Outlier detection

In subsection 3.1, we have cleaned away the abnormal data which deviate too
much from the dataset. However, there are still some outliers. An outlier is a
point which differs significantly from the other points in a subdataset
measured by the same device in a short time. We use the spherical distance
computed by the haversine formula S1984  to identify the outlier. The
haversine formula is illustrated by Fig. 2 and calculates the spherical
distance between the two points $A(lon_{a},\,lat_{a})$ and
$B(lon_{b},\,lat_{b})$ with the coordinate $(longitude,\,latitude)$ as
follows:

$\displaystyle d_{AD}$ $\displaystyle=2Rsin(\Delta lon/2)cos(lat_{a}),$
$\displaystyle d_{CB}$ $\displaystyle=2Rsin(\Delta lon/2)cos(lat_{b}),$ (2)
$\displaystyle d_{AB}$ $\displaystyle=2R\left|sin^{2}\left(\frac{\Delta
lat}{2}\right)+cos(lat_{a})cos(lat_{b})sin^{2}\left(\frac{\Delta
lon}{2}\right)\right|^{\frac{1}{2}},$

where $\Delta lon=lon_{b}-lon_{a}$, $\Delta lat=lat_{b}-lat_{a}$, and $R$ is
the radius of the Earth.

Figure 2: The diagram of two points in a three-dimensional space.

Consequently, we estimate the diameter of a subdataset as follows:

$\displaystyle d_{max}=\bar{v}\times t,$ (3)

where $\bar{v}$ is the mean velocity, and $t$ is the total measuring time of
the subdataset. On the other hand, each point has a distance vector with other
points. If over $50\%$ elements of the distance vector are greater than
$d_{max}$, we regard this point as an outlier.

###  3.3 Data correction

In this subsection, we describe the procedure of data correction and it is
also the key step of our positioning method. This step is to predict the
relatively accurate altitudes of the outliers. As mentioned in section 3.2,
the data distribution is roughly similar when the data are measured by the
same device. Under this assumption, the altitudes of the subdataset are
classified into different classes (labels). Thus, we encounter the multi-class
classification problem.

####  3.3.1 The multi-class classification problem

The outliers of the subdataset have been found with the method in section 3.2.
Thus, we select the data except outliers as a training dataset. The input
training dataset is composed of $N$ pairwise points
$(X_{n},\,h_{n})\,(n=1,\,2,\,\dots,N)$, where $X_{n}$ is the feature vector of
the $n$-th point and $h_{n}$ is the corresponding altitude. Denote $h_{min}$
and $h_{max}$ as the minimum altitude and the maximum altitude, respectively.
Parameter $\delta\,(h_{min}<\delta<h_{max}$) is the quantization step of
altitude. Then, for a given altitude $h$, its corresponding class $k$ is
computed as follows:

$\displaystyle k=\left\lceil{\frac{h-h_{min}}{\delta}}\right\rceil+1,$

where $h_{min}\leq h\leq h_{max}$,$\left\lceil{\cdot}\right\rceil$ is a
function which will round the value toward positive infinity. When the
predicted class of a point is obtained, we take the average altitude of its
corresponding interval as the predicted altitude and which is computed by the
following formula:

$\displaystyle
h^{p}_{k}=\left(k-\frac{1}{2}\right)\delta+h_{min},\;k=1,\,2,\,\ldots,\,K,$
(4)

Thus, after the above transformation procedure, the data correction problem is
converted into a multi-class classification problem (see Table 1, where $K$
represents the number of classes and
$K=\lceil(h_{max}-h_{min})/\delta\rceil$+1).

Table 1: The class, its corresponding interval and predicted altitude. Class | Interval (meter) | Predicted altitude (meter)
---|---|---
1 | $h_{min}\sim\delta+h_{min}$ | $\frac{1}{2}\delta+h_{min}$
2 | $\delta+h_{min}\sim 2\delta+h_{min}$ | $\frac{3}{2}\delta+h_{min}$
⋮ | ⋮ | ⋮
k | $(k-1)\delta+h_{min}\sim k\delta+h_{min}$ | $\left(k-\frac{1}{2}\right)\delta+h_{min}$
⋮ | ⋮ | ⋮
K | $(K-1)\delta+h_{min}\sim K\delta+h_{min}$ | $\left(K-\frac{1}{2}\right)\delta+h_{min}$

####  3.3.2 The weighted multinomial logistic regression model

Logistic Regression (LR) is a machine learning method and widely used to the
binary classification problem C2006 . The MLR method extends the binary LR
method to the multiple classification problem. For the MLR model, each class
has its parameter vector. According to the parameter vector and the data
feature vector, the MLR method determines the classification of the data. In
the positioning application scenario, every feature vector consists of time,
longitude, latitude, air pressure value and speed.

The training process of the MLR model needs to obtain the parameter
$\omega_{k}$ of the $k$-th class via solving the the maximum likelihood
function ZLC2012 , where $k=1,\,2,\,\cdots,K$. The conditional probability of
the feature vector $X$ belonging to the class $Y$ is given by the following
formula:

$\displaystyle
P(Y=k|X=x)=\frac{e^{\omega_{k}^{T}x}}{\sum^{K}_{i=1}e^{\omega_{i}^{T}x}},\;k=1,2,\cdots,K.$
(5)

Then, the MLR method predicts the data category $k^{\ast}$ via solving the
following maximum problem:

$\displaystyle
k^{\ast}\in\mathop{argmax}\limits_{k\in\\{1,\,2,\,\ldots,\,K\\}}P(Y=k|X=x).$
(6)

After the data preprocessing of the previous steps, we obtain the training
dataset, which consists of $N$ pairwise points
$(X_{n},\,Y_{n})\,(n=1,\,2,\,\ldots,\,N)$, where $X_{n}$ represents the data
feature vector and $Y_{n}$ represents its corresponding data class. According
to formula (5) and the independent assumption of the multivariate
distribution, we obtain the likelihood function as follows:

$\displaystyle\prod^{N}_{n=1}P\left(Y=Y_{n}|X=X_{n}\right)=\prod^{N}_{n=1}\left(\frac{e^{\omega_{Y_{n}}^{T}X_{n}}}{\sum^{K}_{k=1}e^{\omega_{k}^{T}X_{n}}}\right).$
(7)

Taking the logarithm of the two sides of formula (7), we obtain the following
log-likelihood function:

$\displaystyle
log\left(\prod^{N}_{n=1}P\left(Y=Y_{n}|X=X_{n}\right)\right)=\sum^{N}_{n=1}\left(\omega_{Y_{n}}^{T}X_{n}-log\left({\sum^{K}_{k=1}e^{\omega_{k}^{T}X_{n}}}\right)\right).$
(8)

Since the value of expression (8) is less than zero, we define function
$f(\Omega)$ as

$\displaystyle
f(\Omega)=\sum^{N}_{n=1}\left(-\omega_{Y_{n}}^{T}X_{n}+log\left({\sum^{K}_{k=1}e^{\omega_{k}^{T}X_{n}}}\right)\right),$
(9)

where $\Omega=[\omega_{1},\,\omega_{2},\,\ldots,\,\omega_{K}]$. Then, we
obtain the maximum likelihood estimation $\Omega^{\ast}$ of parameter matrix
$\Omega$ via solving the following optimization problem:

$\displaystyle\Omega^{\ast}=\mathop{argmin}\limits_{\Omega}\ f(\Omega).$ (10)

Since the training dataset is separable, the value of function $f(\Omega)$ can
be made arbitrarily close to zero via multiplying $\Omega$ by a large value
KS2016 . In order to maintain the finiteness of $\Omega$, we obtain the
parameter matrix $\Omega^{\ast}$ by solving its regularized problem of problem
(9) as follows:

$\displaystyle\Omega^{*}=\mathop{argmin}\limits_{\Omega}\;\left(f(\Omega)+\lambda\eta(\Omega)\right),$
(11)

where $\lambda>0$ is the regularized parameter and the regularized function
$\eta(\Omega)$ is convex and non-smooth. For this convex optimization problem,
there are many efficient optimization methods to tackle it such as the quasi-
Newton BFGS method (p. 198, NW1999 ). Once the MLR model has been trained, we
can predict the data category via solving the maximum problem (6).

We denote $\mathrm{I}=\\{1,\,2,\,\ldots,\,I\\}$ as the index set of the
feature vector $X$, where $I$ represents the dimension of the feature vector
$X$. Select randomly $r$ features from $I$ features and record the index of
selected features as the subset $\mathrm{S}$ of the index set $\mathrm{I}$.
Since the $\ell_{1}$ regularizer is easier to obtain a sparse solution than
the $\ell_{2}$ regularizer, we define a group-$\ell_{1}$-regularizer as

$\displaystyle\eta_{\mathrm{S}}(\Omega)=\sum_{i\in\mathrm{S}}\|[\Omega]_{\mathrm{I}_{i}}\|_{1},$
(12)

where $[\Omega]_{\mathrm{I}_{i}}$ is the $\mathrm{I}_{i}$-th row of parameter
matrix $\Omega$, and $\|x\|_{1}=\sum_{i=1}^{m}|x_{i}|\;\text{for
vector}\;x\in\Re^{m}$. Thus, the problem (11) is written as the following
group-sparse problem:

$\displaystyle\min_{\Omega}\;\left(f(\Omega)+\lambda\eta_{\mathrm{S}}(\Omega)\right).$
(13)

If the parameter $\lambda$ is suitably selected, the solution $\Omega^{\ast}$
of problem (13) will be group-row-sparse KCFH2005 .

After $L$ operations as the procedure above, we obtain $L$ parameter matrices
$\Omega^{\ast}_{1}$, $\Omega^{\ast}_{2},\,\ldots,\,\Omega^{\ast}_{L}$.
Multiply the $L$ parameter matrices $\Omega^{\ast}_{l}\,(l=1,\,2,\,\ldots,L)$
by their corresponding sub-features, then we obtain the predicted categories
$k_{l}^{\ast}\,(l=1,\,2,\,\ldots,\,L)$ with formulas (5)-(6) and its predicted
altitudes $h_{l}^{p}\,(l=1,\,2,\,\ldots,L)$ with formula (4) as follows:

$\displaystyle
k^{\ast}_{l}=\mathop{argmax}\limits_{k\in\\{1,\,2,\,\ldots,\,K\\}}\left(\omega_{k}^{\ast
l}\right)^{T}[X]_{\mathrm{S}_{l}},\;\text{and}\;h^{p}_{l}=\left(k^{\ast}_{l}-\frac{1}{2}\right)\delta+h_{min},\;l=1,\,2,\,\ldots,\,L,$
(14)

where $[X]_{\mathrm{S}_{l}}$ represents the sub-features selected from the
feature vector $X$ and the $i$-th element of $[X]_{\mathrm{S}_{l}}$ equals
$X(\mathrm{S}_{l}(i))$, $\omega_{k}^{\ast l}$ is the $k$-th element of matrix
$\Omega_{l}$.

Compute $L$ absolute errors between the original altitude $h$ and the $l$-th
predicted altitude $h_{l}^{p}(l=1,\,2,\,\ldots,\,L)$ as follows:

$\displaystyle Err_{l}=\left|h-h_{l}^{p}\right|,\;l=1,\,2,\,\ldots,L.$ (15)

Then, we obtain the weighted predicted altitude of the feature vector as
follows:

$\displaystyle h^{\ast}=\sum^{L}_{l=1}w_{l}h_{l}^{p},$ (16)

where the weighted coefficients $w_{l}\,(l=1,\,2,\,\ldots,\,L)$ are computed
by the following formula:

$\displaystyle
w_{l}=\frac{Err_{l}}{\sum^{L}_{l=1}Err_{l}},\;l=1,\,2,\,\ldots,\,L.$ (17)

According to the above discussions, we give the weighted multinomial logistic
regression method for the vertical position problem in Algorithm 1.

Algorithm 1 The WMLR method for the vertical positioning problem

0: the training data $(X_{n},h_{n}),n=1,\,2,\,\ldots,\,N$; the test feature
vector $X$ and its corresponding altitude $h$.

0: the predicted altitude $h^{\ast}$ of the feature vector $X$.

1: Given the regularized parameter $\lambda$, the dimension $r$ of the sub-
feature vector, the quantization step $\delta$ of altitude, the number of the
group-sparse operations $L$.

2: for $l=1,\,2,\,\ldots,\,L$ do

3: Select randomly $r$ features from every feature vector of the training
dataset and denote its corresponding index set of $r$ features as
$\mathrm{S}_{l}$.

4: Obtain the $l$-th regression coefficient matrix $\Omega^{\ast}_{l}$ via
solving the optimization problem
$\Omega_{l}^{\ast}=\mathop{argmin}\limits_{\Omega}\;\left(f(\Omega)+\lambda\eta_{\mathrm{S}_{l}}(\Omega)\right)$,
where $f(\Omega)$ is defined by equation (11) and
$\eta_{\mathrm{S}_{l}}(\Omega)$ is defined by equation (12).

5: Obtain the predicted category $k^{\ast}_{l}$ and the $l$-th predicted
altitude $h_{l}^{p}$ of the feature vector $X$ via solving problem (14).

6: Compute the absolute error $Err_{l}$ between the original altitude and the
predicted altitude of the feature vector $X$ from equation (15).

7: end for

8: Compute $L$ weighted coefficients $w_{l}\,(l=1,\,2,\,\ldots,\,L)$ from
equation (17).

9: Obtain the weighted predicted altitude $h^{\ast}$ of the feature vector $X$
from equation (16).

##  4 Numerical experiments

In this section, we compare the MLR method, the WMLR method (Algorithm 1) and
the SVM method (coded by C. Chang and C. Lin, CL2013 ) for the vertical
positioning problem. The programs are performed under the MATLAB environment
MATLAB .

The raw dataset is provided by Huawei Technologies Company and collected by
different user equipments. From Fig. 3, we find that there are 12796 UserIds
and the number of data collected by each UserId is different. In the dataset,
each piece of data includes time, longitude, latitude, speed, altitude and
some data also contain barometric pressure value. The measurement time of the
experiment dataset spans almost three months from October 5 to December 25,
2018. The air pressure is relatively high because the temperature is
relatively low in that season. Except for null values, the data type is
numeric.

Figure 3: The data volume of the corresponding userID of the raw dataset

Since the raw dataset contains many null and abnormal values, we exclude those
null and abnormal values with the method in subsection 3.1. Table 2 presents
the statistical results of the cleaned data. From Table 2, we find that the
distribution of data is not Gaussian. Thus, we standardize and normalize the
data. After the data cleaning and normalization, we obtain a training set,
every data element of which includes time, speed, longitude, latitude,
pressure. We divide the dataset into two parts, i.e. $70\%$ data for training
and $30\%$ data for testing.

Then, in order to verify the effect of Algorithm 1 (the WMLR method), we
compare the performance of the MLR method, Algorithm 1, and the SVM method for
the cleaned data. For Algorithm 1, we set the regularized parameter
$\lambda=10^{-3}$, the quantization step $\delta=4$, the length of the group-
sparse feature $r=4$ and $L=C_{5}^{4}=5$. The numerical results are put in
Table 3 and Fig. 4. Table 3 is the statistical results of the vertical
positioning accuracy predicted by three methods. From Table 3, we find that
the vertical positioning accuracy is improved from 5.9 meters (the MLR
method), 5.4 meters (the SVM method) to 5 meters (the WMLR method) for $67\%$
test points. Fig. 4 is the cumulative distribution function of the positioning
accuracy. From Fig. 4, we find that the positioning error of WMLR is less than
that of the SVM method and the WLR method when the cumulative probability is
less than 90%, and the positioning accuracy of the SVM method is the best when
the cumulative probability is greater than 90%.

Table 2: The statistical results of the cleaned data. | longitude | latitude | speed | pressure | label | altitude
---|---|---|---|---|---|---
mean | 121.5767 | 31.2595 | 5.8808 | 1021.3788 | 0.9181 | 22.9314
std | 0.0030 | 0.0020 | 6.7051 | 1.2559 | 0.2742 | 10.9594
min | 121.5708 | 31.2566 | 0.0000 | 1017.1787 | 0.0000 | 0.0534
25% | 121.5742 | 31.2579 | 1.0000 | 1020.5680 | 1.0000 | 15.7657
50% | 121.5765 | 31.2590 | 3.0000 | 1021.3281 | 1.0000 | 20.1303
75% | 121.5792 | 31.2610 | 10.0000 | 1022.3744 | 1.0000 | 28.5893
max | 121.5820 | 31.2653 | 26.0000 | 1024.0759 | 1.0000 | 78.1991
Table 3: Vertical positioning accuracies (m) of the MLR, WMLR and SVM methods. | Min | Max | Mean | Median | Std | $67\%$ | $90\%$
---|---|---|---|---|---|---|---
MLR | 0.0211 | 48.8268 | 5.9795 | 4.4133 | 6.7941 | 5.9705 | 11.7054
WMLR | 0.0211 | 31.9072 | 4.6628 | 3.2539 | 3.2539 | 5.0216 | 10.1085
SVM | 0.0211 | 25.2855 | 4.9508 | 3.9297 | 4.0743 | 5.4383 | 10.3968
Figure 4: The comparison of three different methods.

##  5 Conclusion and future works

In this paper, a vertical positioning method with GPS information and the air
pressure values is proposed. Firstly, we clean the missing and abnormal data.
Then, according to the spherical distance matrix between points, we identify
and exclude outliers. Consequently, we divide the cleaned data into two parts,
i.e. $70\%$ data for training and $30\%$ data for testing. Based on the
cleaned data, we compare the performances of the MLR method, the WMLR method
(Algorithm 1), and the SVM method for this vertical positioning problem. The
numerical results show that the vertical positioning accuracy is improved from
5.9 meters (the MLR method), 5.4 meters (the SVM method) to 5 meters (the WMLR
method). Therefore, the WMLR method has some improvements of the positioning
accuracy for this vertical positioning problem.

The appealing positioning technology based on the WMLR method is that this
method does not rely on the empirical pressure-height formula and it can
automatically adjust the parameter matrix according to the local area. The
integration of the MLR method and the weighted technique considers all
training points such that it smoothes the noise to get a better prediction.
For the WMLR method, since it exists the quantization step, it will result in
enlarging the positioning error when the point is the misclassification, which
is a problem to be solved in the future work. Besides, due to the
heterogeneity of user equipments and the complexity of the real environment,
there are some room of improvement on the vertical positioning accuracy of the
WMLR method.

## Financial and Ethical Disclosures

  * $\bullet$

Funding: This work was supported in part by Grant 61876199 from National
Natural Science Foundation of China, Grant YBWL2011085 from Huawei
Technologies Co., Ltd., and Grant YJCB2011003HI from the Innovation Research
Program of Huawei Technologies Co., Ltd..

  * $\bullet$

Conflict of Interest: The authors declare that they have no conflict of
interest.

## References

  * (1) A. B. Adege, H. Lin, G. B. Tarekegn, Y. Y. Munaye and L. Yen, _An indoor and outdoor positioning using a hybrid of support vector machine and deep neural network algorithms_ , Journal of Sensors, 2018, 1-12 (2018).
  * (2) S. Alaee, A. Abdoli, C. Shelton, A. C. Murillo, A. C. Gerry and E. Keogh, _Features or shape? Tackling the false dichotomy of time series classification_ , arXiv preprint, https://arxiv.org/abs/1912.09614, (2019).
  * (3) N. Ali, B. Zafar, M. K. Iqbal, M. Sajid, M. Y. Younis, S. H. Dar, M. T. Mahmood and I. H. Lee, _Modeling global geometric spatial information for rotation invariant classification of satellite images_ , PLoS ONE, 14 (7), 1-24 (2019).
  * (4) S. Burgess, K. Åström, M. Högström and B. Lindquist, _Smartphone positioning in multi-floor environments without calibration or added infrastructure_ , 2016 International Conference on Indoor Positioning and Indoor Navigation (IPIN), IEEE (2016).
  * (5) C. Chang and C. Lin, _Training $\nu$-support vector classifiers: theory and algorithms_, Neural Computation, 13, 2119-2147 (2001).
  * (6) C. Chang and C. Lin, _LIBSVM: a library for support vector machines_ , the software package available at https://www.csie.ntu.edu.tw/~cjlin/libsvm/, 2013.
  * (7) B. Christopher M, _Pattern Recognition and Machine Learning_ , Springer, New York, USA, 2006.
  * (8) A. Chriki, H. Touati and H. Snoussi, _SVM-based indoor localization in wireless sensor networks_ , 2017 13th International Wireless Communications and Mobile Computing Conference, 1144-1149 (2017).
  * (9) H. Du, C. Zhang, Q. Ye, W. Xu, P. L. Kibenge and K. Yao, _A hybrid outdoor localization scheme with high-position accuracy and low-power consumption_ , EURASIP Journal on Wireless Communications and Networking, 2018 (4), 1-13 (2018).
  * (10) P. Gupta, S. Bharadwaj, S. Ramakrishnan and J. Balakrishnan, _Robust floor determination for indoor positioning_ , 2014 Twentieth National Conference on Communications (NCC), Kanpur, 1-6 (2014).
  * (11) T.-Y. He, X.-L. Luo and Z.-H. Liu, _A probabilistic indoor localization algorithm based on restricted Boltzmann machine_ , Proceedings of 2017 IEEE 2nd Advanced Information Technology, Electronic and Automation Control Conference, 1364-1368 (2017).
  * (12) J. Huang and C. Tsai, _Improve GPS positioning accuracy with context awareness_ , 2008 First IEEE International Conference on Ubi-Media Computing, Lanzhou, 94-99 (2008).
  * (13) M. Islam and J. Kim, _An effective approach to improving low-cost GPS positioning accuracy in real-time navigation_ , The Scientific World Journal, 2014, 1-8 (2014), http://dx.doi.org/10.1155/2014/671494.
  * (14) K. Kayabol, _Approximate sparse multinomial logistic regression for classification_ , IEEE Transactions on Pattern Analysis and Machine Intelligence, 42 (2), 490-493 (2019).
  * (15) T. Kim and S. Wright, _PMU placement for line outage identification via multiclass logistic regression_ , IEEE Transactions on Smart Grid, 9 (1), 122-131 (2016).
  * (16) B. Krishnapuram, L. Carin, M. A. T. Figueiredo and A. J. Hartemink, _Sparse multinomial logistic regression: Fast algorithms and generalization bounds_ , IEEE Transactions on Pattern Analysis and Machine Intelligence, 27 (6), 957-968 (2005).
  * (17) Z.-H. Liu, X.-L. Luo and T.-Y. He, _Indoor positioning system based on the improved W-KNN algorithm_ , Proceedings of 2017 IEEE 2nd Advanced Information Technology, Electronic and Automation Control Conference, 1355-1359 (2017).
  * (18) L. Meier, S. V. D. Geer and P. Bühlmann, _The group Lasso for logistic regression_ , Journal of the Royal Statistical Society, 70 (1), 53-71 (2008).
  * (19) MATLAB 9.6.0 (R2019a), The MathWorks Inc., http://www.mathworks.com, 2019.
  * (20) J. Nocedal and S. J. Wright, _Numerical Optimization_ , Springer-Verlag, 1999.
  * (21) C. Ordóñez, J. R. Rodríguez-Pérez, J. J. Moreira, J. M. Matías and E. Sanz-Ablanedo, _Machine learning techniques applied to the assessment of GPS accuracy under the forest canopy_ , Journal of Surveying Engineering, 137, 140-149 (2011).
  * (22) M. Rohani, D. Gingras and D. Gruyer, _A novel approach for improved vehicular positioning using cooperative map matching and dynamic base station DGPS concept_ , IEEE Transactions on Intelligent Transportation Systems, 17 (1), 230-239 (2015).
  * (23) R. W. Sinnott, _Virtues of the haversine_ , Sky and Telescope, 68 (2), 158-159 (1984).
  * (24) Y. H. Wang, H. Li, X.-L. Luo, Q. M. Sun, J. N. Liu, _A 3D fingerprinting positioning method based on cellular networks_ , International Journal of Distributed Sensor Networks, 1-9 (2014), http://dx.doi.org/10.1155/2014/248981.
  * (25) P. Williams, _Interactive Statistics for the Behavioral Sciences_ , Sinauer Associates, Inc. Publishers, Sunderland, Massachusetts USA, 2004.
  * (26) H. Xia, X. Wang, Y. Qiao, J. Jian and Y. Chang, _Using multiple barometers to detect the floor location of smart phones with built-in barometric sensors for indoor positioning_ , Sensors, 15(4), 7857-7877 (2015).
  * (27) V. Zaliva and F. Franchetti, _Barometric and GPS altitude sensor fusion_ , 2014 IEEE International Conference on Acoustics, Speech and Signal Processing, 1-5 (2014).
  * (28) J. Zhu, X.-L. Luo and D. Chen, _Maximum likelihood scheme for fingerprinting positioning in LTE system_ , 2012 IEEE 14th International Conference on Communication Technology, 428-432 (2012).
  * (29) H. Zou, B. Huang, X. Lu, H. Jiang and L. Xie, _A robust indoor positioning system based on the procrustes analysis and weighted extreme learning machine_ , IEEE Transactions on Wireless Communications, 15 (2), 1252-1266 (2016).


# Deep Learning-based Group Causal Inference in Multivariate Time-series

Wasim Ahmad, Maha Shadaydeh, Joachim Denzler

###### Abstract

Causal inference in a nonlinear system of multivariate time series is
instrumental in disentangling the intricate web of relationships among
variables, enabling us to make more accurate predictions and gain deeper
insights into real-world complex systems. Causality methods typically identify
the causal structure of a multivariate system by considering the cause-effect
relationship of each pair of variables while ignoring the collective effect of
a group of variables or interactions involving more than two-time series
variables. In this work, we test model invariance by group-level interventions
on the trained deep networks to infer causal direction in groups of variables,
such as climate and ecosystem, brain networks, etc. Extensive testing with
synthetic and real-world time series data shows a significant improvement of
our method over other applied group causality methods and provides us insights
into real-world time series. The code for our method can be found at:
https://github.com/wasimahmadpk/gCause.

## Introduction

Group-based causal inference investigates causal relationships within specific
groups of individuals, entities, or units. It is particularly relevant when
studying complex systems with interconnected components such as in climate
(Molotoks et al. 2020) and brain networks (Faes et al. 2022), allowing
researchers to investigate how variables within distinct groups contribute to
observed outcomes or behaviors. In this paper, we explore the fundamental
concept behind group-based causal inference in multivariate time series and
its growing significance in diverse disciplines. To this end, we present our
approach to testing the causal link in groups of time series, see Figure 1,
which provides a comprehensive understanding of how these groups interact. The
proposed method allows the testing of bi-directional causal links. Our method
builds on our previous work (Ahmad, Shadaydeh, and Denzler 2022), which
exploits the model invariance property through Knockoffs (Barber and Candès
2015, 2019; Barber, Candès, and Samworth 2020) interventions for pairs of
variables in deep networks for causal estimation. However, here we emphasises
the identification of causal interaction in groups of variables and is the
first group causality method to the best of our knowledge which utilize deep
learning to learn complex nonlinear relation. The model invariance property
refers to the invariant behavior of the model in different settings when its
causal predictors are observed (Peters, Bühlmann, and Meinshausen 2016). We
use DeepAR (Salinas et al. 2020) to model multivariate time series which has
the potential to learn complex interactions in nonlinear data. Since deep
networks cannot handle missing variables or adjust to out-of-distribution
data, we use Knockoff variables for interventions on the trained deep
networks. Knockoffs are in-distribution, uncorrelated copies of the original
data with similar covariance structure.

Figure 1: Interactions among groups of multivariate time series data of
various dimensions.

To demonstrate the robustness and practical utility of our approach, we
evaluate its performance on synthetically generated time series and real-world
datasets. We compare our method with other group causality methods, i.e.,
Vanilla-PC (Janzing, Hoyer, and Schölkopf 2009), Trace method (Zscheischler,
Janzing, and Zhang 2012) and 2GVecCI (Wahl, Ninad, and Runge 2023). As a real-
world application, we mainly focus on the climate system and assessing the
direction of connections in brain networks. This involves considering the
internal processes over a period of time involving interactions within and
between different climate subsystems, i.e., ocean, land, ecosystem,
atmosphere, etc. (Latif and Keenlyside 2009). For example, how climate and
ecosystem interact (Malhi et al. 2020; Sefidmazgi and Sefidmazgi 2020; Korell
et al. 2020), or how major climate phenomena such as the El Niño Southern
Oscillation (ENSO) affect different global regions (Nowack et al. 2020). These
highly complex interactions provide the basis for understanding the profound
impacts of environmental change. We conducted experiments on the FLUXNET data
(Pastorello et al. 2020), which contains measurements from different sites
that capture climate-ecosystem dynamics, where we aim to uncover the causal
pathways between climate and ecosystem. We have also tested our method on a
climate science dataset (ENSO 3.4), where we aim to estimate the impact of
temperature patterns in the tropical Pacific on British Columbia. Moreover, we
apply our method to simulated fMRI time series to identify connections in the
brain network and indicate their directionality.

## Related Work

There are many methods for time series causal inference with applications in
various fields that rely on certain assumptions that limit their applicability
(Yao et al. 2020). These approaches treat each system variable separately and
estimate its causal relation to other variables, overlooking the group or
collective influence, e.g., causal relation in brain regions (Siddiqi et al.
2022), ecosystem and climate subsystems, etc. The work of (Besserve et al.
2018) presents a group theoretical framework to assess the relationship
between cause and mechanism using group transformations. Vanilla-PC (Janzing,
Hoyer, and Schölkopf 2009) estimates the direction of a causal link in two
groups of time series; however, the method suffers from high dimensionality in
terms of accuracy and computation time because it runs a higher number of
conditional independence tests. The Trace method is introduced by
(Zscheischler, Janzing, and Zhang 2012) to infer causal direction in two
linearly interacting groups which works faster, however, its performance
degrades in case of high nonlinear relationships in groups. The authors of
(Wahl, Ninad, and Runge 2023) present 2GVecCI for uni-directional causality
between two groups where they combine the constraint-based approach for causal
discovery with sparsity measures of the internal causal structure of the
groups. The 2GVecCI method is based on two principles for causal relations
between groups of time series. First, conditioning on the cause group does not
lead to new conditional dependencies within the effect group. Second,
conditioning on the effect group does not delete conditional dependencies
within the cause group.

The authors of (Faes et al. 2022) assess higher-order interactions in networks
of processes through time and frequency domain analysis. The work of (Sato et
al. 2010) uses Granger causality for a set of time series to analyze the
connectivity between brain regions, assuming that sets are linearly related to
each other while real-world data are highly nonlinear. Approaches that depend
on the non-Gaussian nature of noise for drawing conclusions about causal
direction, for example, LinGaM (Shimizu et al. 2006), become less effective
due to the influence of averaging caused by the central limit theorem, which
tends to make the averaged noise distribution more Gaussian. Furthermore,
methods for group causality aim to transform data by aggregation or
dimensionality reduction, where a significant amount of information may be
lost and, more importantly, the existing causal structure (Spirtes, Glymour,
and Scheines 2000) is disrupted. However, our work utilizing deep networks has
the potential to learn complex relationships in high-dimensional data.

Figure 2: Schematic diagram of group causal discovery method where we model
the complex relation among all variables in the system with deep learning
approximator. For group-level causal inference, we implement model invariance
testing by applying intervention to a group of interest and estimate its
influence on the target group.

For brain data, the detection of networks and their directions using fMRI time
series is of high interest. The work of (Smith et al. 2011) presents various
correlation and causality-based method to identify the association between
brain network nodes and assign a direction to them. They generate realistic
simulated fMRI data for a wide range of underlying networks and experimental
protocols to compare different connectivity estimation approaches.

Figure 3: Distribution shift of a. ecosystem group
$G_{\text{E}}=\\{GPP,R_{\text{eco}}\\}$ in response to intervention on climate
group $G_{\text{C}}=\\{T,R_{g}\\}$. b. British Columbia temperature in
response to intervention on ENSO temperature data. c. Brain network $N_{2}$ in
fMRI time series after group-level intervention on network $N_{1}$.

## Method

In this section, we present the formulation of our method for group causal
discovery assessing the causal relationships within sets of $N$-variate system
of time series denoted as $Z^{N}=\\{Z_{1},\dots,Z_{N}\\}$. These variables are
organized into $G$ groups or subsets of variables, represented as
$X^{G}=\\{X_{1},\dots,X_{G}\\}$, which can be considered as indicative of the
behavior of a network comprising $G$ distinct dynamic subsystems. Referring to
the real-world applications, the groups $X_{i},i=1,\dots,G$ may represent
subsets of the climate system, brain networks, etc., where each group process
$Z_{t,t=1,\dots,n}$, may represent the intra-atmospheric interaction or neural
activity of a brain region. Each of these groups $X_{i}:Z\subset Z^{N}$ has a
certain dimension $G_{i}$, such that the total number of variables $N$ can be
expressed as $N=\sum_{i=1}^{G}G_{i}$. The method allows the assessment of uni-
and bi-directional causal links in any number of groups. For two causally
related groups $X_{i},X_{j},i,j=1,\dots,G$ where $i\neq j$, the possible links
can be $X_{i}\overset{\text{causes}}{\longrightarrow}X_{j}$, or
$X_{j}\overset{\text{causes}}{\longrightarrow}X_{i}$ or
$X_{i}\overset{\text{causes}}{\longleftrightarrow}X_{j}$.

### Assumptions

We infer causal relationships in groups of nonlinear multivariate time series
by making the following assumptions:

  * -

Stationarity: The variables $Z^{N}$ represent a stochastic stationary process.

  * -

Causal Sufficiency: The set of observed variables
$Z^{N}=\\{Z_{1},\dots,Z_{N}\\}$ contain all common causes in $Z^{N}$, i.e., no
hidden confounders.

  * -

Model Invariance: The causal structure in the stochastic processes $Z^{N}$
remains consistent across different interventional enviroments.

The presented method for group causality in time series builds on our previous
work (Ahmad, Shadaydeh, and Denzler 2022) where we model the complex
interactions in $N$-variate time series by training deep networks and apply
model invariance testing through group-level interventions for inferring
causal direction in groups. A model is invariant if, in the presence of its
causal predictors, the distribution of its output residuals does not change
across interventional environments (Peters, Bühlmann, and Meinshausen 2016).
The schematic diagram for our method is shown in Figure 2, where we apply
invariance testing to deep networks trained on sets of time series in order to
infer causality in groups. Through deep autoregressive models (Salinas et al.
2020) which takes $Z^{N}=\\{Z_{1},\dots,Z_{N}\\}$ as input, we model the
conditional distribution
$P(Z_{i,t_{0}:T}|Z_{i,1:t_{0}-1},Z_{-i,1:t_{0}-1}),i=1,\dots,N$ of the future
of all system variables given their past values. Here $Z_{i}$ represents node
of the target group and $Z_{-i}$ refers to the nodes in the rest of the groups
in $X^{G}$. It utilizes recurrent neural networks (RNNs) to generate
probabilistic output in terms of $\mu$ and $\sigma^{2}$ at each time point
where the mean represents the central estimate, and the variance represents
the uncertainty of the system. For all variables, we obtain distribution $R$
of their residuals $e^{n}=e_{1},\dots,e_{n}$ for $n$ forecast windows where
$e=\frac{1}{T}\sum_{t=1}^{T}|\frac{Z_{t}-Z_{pred}}{Z_{t}}|$ and $T$ represents
forecast horizon. Pertaining to the invariance property, the model response
$Z_{i}$, which is the target group of variables, does not change in different
settings as long as we don’t perturb the cause group.

### Knockoff Intervention

The effectiveness of our previously proposed Knockoff intervention (Barber and
Candès 2015, 2019; Barber, Candès, and Samworth 2020) compared to other
intervention methods is demonstrated on node-to-node causal relationship in
multivariate time series (Ahmad, Shadaydeh, and Denzler 2021, 2022). Let’s say
$X_{j}$ represent the group of output variables of the trained deep network,
and $X_{i}$ denotes the input group where $i\neq j$. The causal effect of
group $X_{i}$ on the network output group $X_{j}$ under intervention is
denoted as:

$E(X_{j}|do(X_{i}=\tilde{X_{i}})$ (1)

Figure 4: a. Illustration of exchangeability property
$(Z_{1},Z_{2})\overset{d}{=}(Z_{1},\tilde{Z}_{2})$ of knockoff variables b.
Shows the correlation matrix $\Sigma_{Z\tilde{Z}}$ of the original variables
$Z$ and knockoffs $\tilde{Z}$ where the sub-matrix for the correlation within
the generated knockoffs $\Sigma_{\tilde{Z}_{1}\tilde{Z}_{2}}$ is similar to
that of correlation in original variables $\Sigma_{Z_{1}Z_{2}}$ , enclosed in
red squares. While the variable-wise correlation is minimized with their
respective knockoff copies, i.e.,
$\sigma_{Z_{1}\tilde{Z}_{1}},\sigma_{Z_{2}\tilde{Z}_{2}}\approx 0$, shown in
white squares.

where $\tilde{X_{i}}$ is a knockoff representation or configuration to which
the group of variables $X_{i}$ is set. Note that group $X_{i}$ contains
internal process or nodes $Z_{t},t=1,\dots,n$. Each node within a group is
replaced with the generated knockoff copies. Applying the do-operator on
trained deep networks involves specifying interventions or changes to the
input group of time series to observe the corresponding effects on the target
group. This implies that, instead of observing the natural variation in
$X_{i}$, we are setting it to knockoff version $\tilde{X_{i}}$ in order to
assess the collective causal impact on the network’s output. We choose to use
multivariate Gaussian model within the Knockoff framework which implements
semidefinite programming (SDP) to estimate knockoff parameters given the mean
$\mu_{Z}$ and covariance matrix $\Sigma_{Z}$ of the original variables $Z$. We
generate knockoff copies $\widetilde{Z}_{1},\dots,\widetilde{Z}_{N}$ as an
interventional representation of $Z_{1},\dots,Z_{N}$, which are used to
replace variables during an intervention. The generated knockoff variables are
in-distribution and uncorrelated with the original variables. Moreover, it
fulfills the exchangeability condition
$(Z_{i},\tilde{Z}_{j})\overset{d}{=}(\tilde{Z}_{i},Z_{j}),i,j=1,\dots,N,i\neq
j$ (Barber and Candès 2019) by maintaining same covariance structure as the
original time series, which means that the underlying joint distribution of
the data does not change by replacing a variable with its knockoff copy. The
properties of knockoffs are demonstrated in Figure 4.

Table 1: Performance of the applied methods on synthetic data for varying interaction densities. The presented values are the ratio of the method output and the total number of tests where the Correct inference represents the true causal direction, Wrong is the opposite of the true causal direction, and No inference is reporting the absence of a causal link. |  | Interaction density |
---|---|---|---
Method | Inference | 0.2 | 0.3 | 0.4 | 0.5 | 0.6 | 0.7 | 0.8 | 0.9 | 1.0
gCDMI | Correct | 1.00 | 0.67 | 1.00 | 1.00 | 1.00 | 0.67 | 0.67 | 1.00 | 1.00
Wrong | 0.00 | 0.33 | 0.00 | 0.00 | 0.00 | 0.33 | 0.33 | 0.00 | 0.00
| No inference | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00
Trace | Correct | 0.67 | 1.00 | 0.67 | 0.67 | 0.67 | 0.67 | 0.67 | 1.00 | 0.67
Wrong | 0.33 | 0.00 | 0.33 | 0.33 | 0.33 | 0.33 | 0.33 | 0.00 | 0.33
| No inference | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00
2GVecCI | Correct | 0.67 | 0.00 | 1.00 | 0.33 | 0.67 | 0.33 | 1.00 | 0.67 | 0.00
Wrong | 0.33 | 0.33 | 0.00 | 0.33 | 0.00 | 0.00 | 0.00 | 0.00 | 0.33
| No inference | 0.00 | 0.67 | 0.00 | 0.33 | 0.33 | 0.67 | 0.00 | 0.33 | 0.67
Vanilla-PC | Correct | 0.67 | 0.67 | 0.67 | 0.33 | 0.67 | 0.67 | 0.33 | 0.33 | 0.00
Wrong | 0.00 | 0.00 | 0.33 | 0.33 | 0.00 | 0.33 | 0.67 | 0.67 | 1.00
| No inference | 0.33 | 0.33 | 0.00 | 0.33 | 0.33 | 0.00 | 0.00 | 0.00 | 0.00

### Group Causal Inference

In order to infer causal direction in two groups of variables, i.e.,
$X_{i}\rightarrow X_{j}$, we apply Kolmogorov–Smirnov (KS) test (Smirnov 1939)

$C_{ij}=\sqrt{\frac{qr}{q+r}}sup|R-\widetilde{R}|$ (2)

to evaluate model invariance property by estimating the maximum absolute
difference in the marginal residual distribution $R$ and $\widetilde{R}$ of
variables in $X_{j}$ before and after group-level intervention on $X_{i}$.
Here $\widetilde{R}$ represents the distribution of the residuals
$\widetilde{e}^{n}=\widetilde{e}_{1},\dots,\widetilde{e}_{n}$, obtained while
estimating the counterfactuals
$P(Z_{i,t_{0}:T}|Z_{i,1:t_{0}-1},\widetilde{Z}_{-i,1:t_{0}-1})$. The KS test
uses a significance level $\alpha$ and other parameters like $q,r$ which
represents the number of samples in marginal residual distributions $R$ and
$\widetilde{R}$. The test statistic $C$ is used to calculate $p$ value for
deciding whether a causal link exists or not. If $p_{ij}>\alpha$, the null
hypothesis $H_{0}$: $X_{i}$ does not cause $X_{j}$ is accepted, which means
that the residual distributions of group of variables before and after
intervention are approximately identical across various interventional
settings, i.e., the group invariance property is fulfilled. In the alternate
case, the null hypothesis is rejected, i.e. $H_{1}$: $X_{i}$ causes $X_{j}$.
We test for both causal links $X_{i}\rightarrow X_{j}$ and $X_{j}\rightarrow
X_{i}$ for all possible values of $i,j$ where $i\neq j$. It is important to
mention that if a single node within the group is influenced by its causal
group, it establishes a causal connection between the groups and it doesn’t
necessarily imply that the entire target group is uniformly affected by its
cause.

Table 2: Performance of different group causality methods in identifying causal direction in climate-ecosystem data for various sites. Here $\checkmark$ indicates the presence of the link, while $\times$ indicates the absence of the link and $\rightarrow$ represents the direction of causal relation between groups. We test a bi-directional causal link only for our method. |  | Methods
---|---|---
Sites | Inference | gCDMI | Trace | V-PC | 2GVCI
DE-Hai | $G_{\text{C}}\rightarrow G_{\text{E}}$ | $\times$ | $\times$ | $\times$ | $\times$
$G_{\text{C}}\leftarrow G_{\text{E}}$ | $\times$ | $\checkmark$ | $\times$ | $\times$
| $G_{\text{C}}\leftrightarrow G_{\text{E}}$ | $\checkmark$ | $-$ | $-$ | $-$
| $G_{\text{C}}\nleftrightarrow G_{\text{E}}$ | $\times$ | $\times$ | $\checkmark$ | $\checkmark$
IT-MBo | $G_{\text{C}}\rightarrow G_{\text{E}}$ | $\checkmark$ | $\checkmark$ | $\checkmark$ | $\times$
$G_{\text{C}}\leftarrow G_{\text{E}}$ | $\times$ | $\times$ | $\times$ | $\times$
| $G_{\text{C}}\leftrightarrow G_{\text{E}}$ | $\times$ | $-$ | $-$ | $-$
| $G_{\text{C}}\nleftrightarrow G_{\text{E}}$ | $\times$ | $\times$ | $\times$ | $\checkmark$
FR-Pue | $G_{\text{C}}\rightarrow G_{\text{E}}$ | $\checkmark$ | $\times$ | $\times$ | $\times$
$G_{\text{C}}\leftarrow G_{\text{E}}$ | $\times$ | $\checkmark$ | $\times$ | $\times$
| $G_{\text{C}}\leftrightarrow G_{\text{E}}$ | $\times$ | $-$ | $-$ | $-$
| $G_{\text{C}}\nleftrightarrow G_{\text{E}}$ | $\times$ | $\times$ | $\checkmark$ | $\checkmark$
US-Ton | $G_{\text{C}}\rightarrow G_{\text{E}}$ | $\checkmark$ | $\times$ | $\times$ | $\times$
$G_{\text{C}}\leftarrow G_{\text{E}}$ | $\times$ | $\checkmark$ | $\times$ | $\checkmark$
| $G_{\text{C}}\leftrightarrow G_{\text{E}}$ | $\times$ | $-$ | $-$ | $-$
| $G_{\text{C}}\nleftrightarrow G_{\text{E}}$ | $\times$ | $\times$ | $\checkmark$ | $\times$

## Experiments

### Synthetic Data

To evaluate the performance of our method, we use synthetic data model
$Z_{t}^{j}=\Sigma_{i}f_{i}(Z_{t-k}^{i})+\eta_{t}^{j}$, $i,j={1,...,N},0<k<t$.
The system variables $Z_{j}$ have auto and cross-functional dependencies with
a time delay of $k$. The data model incorporates linear and nonlinear
dependencies $f$, i.e., exponential, polynomial with varying edge densities,
and adds uncorrelated, normally distributed noise $\eta_{t}^{j}$. We
incorporate inter- and intra-group causal links in the generated causal
graphs. The decision of the methods for each edge in the groups is either
correct, wrong or no inference. Our method (gCDMI) achieves better
identification of the correct causal links with fewer wrong detections in
subsystems for all edge densities, i.e., sparse to dense causal graphs, as
compared to other group causality methods as shown in Table 1. The given
values are the ratio of the method decision and the total number of
experiments. It can be noticed that our method reported only fewer wrong
inferences and zero no inference along with the Trace method for all
experiments, while other methods reported no inference in a number of tests
because they are uni-directional and sometimes it is hard to distinguish cause
from effect when the interaction in both directions is not significantly
different. We also found that increasing edge density does not have a clear
impact on the outcome except for Vanilla-PC, where the performance is
degraded. We used multiple group dimensions for each edge density in our
experiments, where we vary the architecture of the deep networks based on
group dimension. All experiments in this paper were conducted on an NVIDIA
GeForce RTX 3060 Ti graphics card, providing reliable performance capabilities
for complex synthetic and real-world data causal analysis. The improved
performance of our method comes at the cost of high computation time due to
its dependence on deep networks to learn complex causal structures in data.

### FLUXNET Data

Here, we aim to perform a causal analysis of environmental time series with
our method by establishing the causal structure in groups or subsystems. We
carry out experiments with FLUXNET2015 dataset (Pastorello et al. 2020), which
is acquired using the eddy covariance technique to measure the cycling of
carbon, water, and energy between the biosphere and atmosphere through
collaborations among many regional networks, with data preparation efforts
happening at site and network levels. For our experiments, we considered
various measurement sites, i.e., Hainich (DE-Hai: Deciduous Broadleaf
Forests), Monte Bondone (IT-MBo: Grasslands), Puechabon (FR-Pue: Evergreen
Broadleaf Forests), and Tonzi Ranch (US-Ton: Woody Savannas) site. The dataset
includes climatic and ecological time series, i.e., global radiation
($R_{\text{g}}$), temperature ($T$), gross primary production ($GPP$),
ecosystem respiration ($R_{\text{eco}}$) for various time scales i.e., half-
hourly, hourly, daily, weekly and so on. We categorized these variables into
two groups: climate group $G_{\text{C}}$ which contains $T$ and $R_{g}$ and
ecosystem group $G_{\text{E}}$, which consists of the ecosystem variables
$GPP$ and $R_{\text{eco}}$. We considered daily sampling which is advantageous
for mitigating the effects of daily patterns that usually undermine the
underlying causal relation in data. The results for all applied causality
methods from various sites are given in Table 2. Our method identified the
presence of a causal link $G_{\text{C}}\rightarrow G_{\text{E}}$ for all sites
except for the DE-Hai site where we obtain a bidirectional link
$G_{\text{C}}\leftrightarrow G_{\text{E}}$ which is an overall better
performance compared to other methods. The bidirectional link at the DE-Hai
site, which is a deciduous broadleaf forest, could be due to climate-ecosystem
strong feedback mechanism or other influential factors that need further
investigation. For DE-Hai, Vanilla-PC and 2GVecCI could not infer any causal
direction while the Trace method detected $G_{\text{E}}\rightarrow
G_{\text{C}}$. For IT-MBo, all methods correctly identified the expected
causal directions except 2GVecCI. As an illustration, we show climatic
influence on the ecosystem for the IT-MBo site from one of our experiments in
Figure 3 (a).

### ENSO Data

Moreover, we performed experiments on the ENSO dataset where we consider
surface temperatures over the ENSO region and British Columbia (BCT) from 1948
to 2021, as a causal effect of temperatures in the tropical Pacific on those
in North America is recognized in climate (Taylor 1998). During an El Niño
event, which is one phase of ENSO, the tropical Pacific Ocean warms up
significantly, disrupt normal weather patterns and influence the climate in
other parts of the world, including the British Columbia region. For
experiments, we adapted data preprocessing, i.e., deseasonalizing, smoothing,
and aggregation from the work of (Runge et al. 2019; Wahl, Ninad, and Runge
2023). We show results for applied group causality methods in Table 3. Both
our method and 2GVecCI identified significant influence of ENSO on BCT for
various grid scales with a fraction of 0.66 correct inferences and 0.34 as no
inferences. While Trace method detected the causal link ENSO $\rightarrow$ BCT
with a fraction of 0.50 correct inferences and 0.50 wrong inferences. Vanilla-
PC could not inferred any causal direction for all grid scales that is
probably because of difficultly in distinguishing data patterns at both
regions. Illustration of the causal influence of ENSO on BCT for one of the
performed tests by our method is given in Figure 3 (b). Which shows the
comparison of the BCT distribution with and without group-level intervention
on the ENSO time series variables. The shift in the counterfactual
distribution for BCT is indicative of the influence of ENSO.

Table 3: Effect of surface temperature at Tropical Pacific ocean (ENSO region) on British Columbia region, computed at various grid scales. The given values represent the decision of the applied methods (in percentage). Inference | gCDMI | Trace | 2GVCI | V-PC
---|---|---|---|---
ENSO $\rightarrow$ BCT | 0.66 | 0.50 | 0.66 | 0.00
ENSO $\leftarrow$ BCT | 0.00 | 0.50 | 0.00 | 0.00
ENSO $\nleftrightarrow$ BCT | 0.34 | 0.00 | 0.34 | 1.00

### fMRI Data

We analyzed simulated fMRI time series data (Smith et al. 2011) to identify
connections in brain network nodes and assign direction to them. These fMRI
time series simulations were based on dynamic causal modelling (DCM) (Friston,
Harrison, and Penny 2003) fMRI forward model. The generated dataset provides
ground truth connectivity graph for variety of network topologies which we use
to evaluate our methods. Here we conducted experiments on a number of
simulated subjects from S2 topology in the dataset which contains 10 nodes
clustered into 2 groups. It has 10 min fMRI sessions for each subject with 3
sec sampling rate, ﬁnal added noise of 1%, and haemodynamic response function
(HRF) variability of ±0.5 which provides time series of 200 data points. The
percentage of methods outcome in terms of correct, wrong and no inference is
given in Table 4. Since our method is data demanding as it relies on deep
networks, the provided size of time series was not long enough for our method
to learn the complex relationship properly. Determining the correct
directionality of connections in fMRI time series is challenging because of
the complex interactions among brain networks. However, our approach inferred
correct direction 56% of the times in brain networks with 21% bidirectional
links which is overall better performance compared to other applied group
causality methods. Results from Trace method were somehow comparable to that
of gCDMI. While 2GVecCI and Vanilla-PC yielded high percentage of no inference
probably because of difficulty in separating cause from effect in a mutual
interaction scenario between groups. The causal influence of brain network
$N_{1}$ on network $N_{2}$ for one subject in S2 topology from simulated fMRI
time series is demonstrated in Figure 3 (c).

Table 4: Identification of the direction of brain network connections using simulated fMRI time series by the applied group causality methods: Ground truth network connection: $N_{1}\rightarrow N_{2}$. The given values represent the decision of the applied methods (in percentage). Inference | gCDMI | Trace | 2GVCI | Vanilla-PC
---|---|---|---|---
$N_{1}\rightarrow N_{2}$ | 0.56 | 0.50 | 0.33 | 0.06
$N_{1}\leftarrow N_{2}$ | 0.17 | 0.39 | 0.17 | 0.16
$N_{1}\leftrightarrow N_{2}$ | 0.21 | - | - | -
$N_{1}\nleftrightarrow N_{2}$ | 0.06 | 0.11 | 0.50 | 0.78

## Conclusion

In this work, we have introduced a deep learning-based method to uncover the
causal interactions in groups or subsets of time series. We model complex
relationships in an $N$-variable system of time series using deep learning
models, and apply invariance testing via group-level intervention to infer
causal direction in groups of variables. Our approach also tests for bi-
directional causal links which signify a mutual influence, suggesting that the
variables have a reciprocal cause-and-effect relationship. While our approach
demonstrated improved performance on both synthetic and real-word time series
data compared to other causality methods, it is worth noting that the high
computational time is a trade-off inherent to its reliance on deep learning.
As future work, we consider estimating the causal interactions in more than
two groups of time series. Moreover, we aim to address the issue of non-
stationarity and hidden confounding (Trifunov, Shadaydeh, and Denzler 2022) to
further improve our method.

## Acknowledgments

This work is funded by the German Research Foundation (DFG) research grant SH
1682/1-1 and the Carl Zeiss Foundation within the scope of the program line
Breakthroughs: Exploring Intelligent Systems for Digitization — explore the
basics (No P2017-01-003), use applications.

## References

  * Ahmad, Shadaydeh, and Denzler (2021) Ahmad, W.; Shadaydeh, M.; and Denzler, J. 2021.  Causal Inference in Non-linear Time-series using Deep Networks and Knockoff Counterfactuals.  In _2021 20th IEEE International Conference on Machine Learning and Applications (ICMLA)_ , 449–454. IEEE.
  * Ahmad, Shadaydeh, and Denzler (2022) Ahmad, W.; Shadaydeh, M.; and Denzler, J. 2022.  Causal Discovery using Model Invariance through Knockoff Interventions.  In _ICML 2022: Workshop on Spurious Correlations, Invariance and Stability_.
  * Barber and Candès (2015) Barber, R. F.; and Candès, E. J. 2015.  Controlling the false discovery rate via knockoffs.  _The Annals of Statistics_ , 43(5): 2055–2085.
  * Barber and Candès (2019) Barber, R. F.; and Candès, E. J. 2019.  A knockoff filter for high-dimensional selective inference.  _The Annals of Statistics_ , 47(5): 2504–2537.
  * Barber, Candès, and Samworth (2020) Barber, R. F.; Candès, E. J.; and Samworth, R. J. 2020.  Robust inference with knockoffs.  _The Annals of Statistics_ , 48(3): 1409–1431.
  * Besserve et al. (2018) Besserve, M.; Shajarisales, N.; Schölkopf, B.; and Janzing, D. 2018.  Group invariance principles for causal generative models.  In _International Conference on Artificial Intelligence and Statistics_ , 557–565. PMLR.
  * Faes et al. (2022) Faes, L.; Mijatovic, G.; Antonacci, Y.; Pernice, R.; Barà, C.; Sparacino, L.; Sammartino, M.; Porta, A.; Marinazzo, D.; and Stramaglia, S. 2022.  A new framework for the time-and frequency-domain assessment of high-order interactions in networks of random processes.  _IEEE Transactions on Signal Processing_ , 70: 5766–5777.
  * Friston, Harrison, and Penny (2003) Friston, K. J.; Harrison, L.; and Penny, W. 2003.  Dynamic causal modelling.  _Neuroimage_ , 19(4): 1273–1302.
  * Janzing, Hoyer, and Schölkopf (2009) Janzing, D.; Hoyer, P. O.; and Schölkopf, B. 2009.  Telling cause from effect based on high-dimensional observations.  _arXiv preprint arXiv:0909.4386_.
  * Korell et al. (2020) Korell, L.; Auge, H.; Chase, J. M.; Harpole, S.; and Knight, T. M. 2020.  We need more realistic climate change experiments for understanding ecosystems of the future.  _Global Change Biology_ , 26(2): 325–327.
  * Latif and Keenlyside (2009) Latif, M.; and Keenlyside, N. S. 2009.  El Niño/Southern Oscillation response to global warming.  _Proceedings of the National Academy of Sciences_ , 106(49): 20578–20583.
  * Malhi et al. (2020) Malhi, Y.; Franklin, J.; Seddon, N.; Solan, M.; Turner, M. G.; Field, C. B.; and Knowlton, N. 2020.  Climate change and ecosystems: Threats, opportunities and solutions.
  * Molotoks et al. (2020) Molotoks, A.; Henry, R.; Stehfest, E.; Doelman, J.; Havlik, P.; Krisztin, T.; Alexander, P.; Dawson, T. P.; and Smith, P. 2020.  Comparing the impact of future cropland expansion on global biodiversity and carbon storage across models and scenarios.  _Philosophical Transactions of the Royal Society B_ , 375(1794): 20190189.
  * Nowack et al. (2020) Nowack, P.; Runge, J.; Eyring, V.; and Haigh, J. D. 2020.  Causal networks for climate model evaluation and constrained projections.  _Nature communications_ , 11(1): 1415.
  * Pastorello et al. (2020) Pastorello, G.; Trotta, C.; Canfora, E.; Chu, H.; Christianson, D.; Cheah, Y.-W.; Poindexter, C.; Chen, J.; Elbashandy, A.; Humphrey, M.; et al. 2020.  The FLUXNET2015 dataset and the ONEFlux processing pipeline for eddy covariance data.  _Scientific data_ , 7(1): 1–27.
  * Peters, Bühlmann, and Meinshausen (2016) Peters, J.; Bühlmann, P.; and Meinshausen, N. 2016.  Causal inference by using invariant prediction: identification and confidence intervals.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 78(5): 947–1012.
  * Runge et al. (2019) Runge, J.; Nowack, P.; Kretschmer, M.; Flaxman, S.; and Sejdinovic, D. 2019.  Detecting and quantifying causal associations in large nonlinear time series datasets.  _Science Advances_ , 5(11): eaau4996.
  * Salinas et al. (2020) Salinas, D.; Flunkert, V.; Gasthaus, J.; and Januschowski, T. 2020.  DeepAR: Probabilistic forecasting with autoregressive recurrent networks.  _International Journal of Forecasting_ , 36(3): 1181–1191.
  * Sato et al. (2010) Sato, J. R.; Fujita, A.; Cardoso, E. F.; Thomaz, C. E.; Brammer, M. J.; and Amaro Jr, E. 2010.  Analyzing the connectivity between regions of interest: an approach based on cluster Granger causality for fMRI data analysis.  _Neuroimage_ , 52(4): 1444–1455.
  * Sefidmazgi and Sefidmazgi (2020) Sefidmazgi, M. G.; and Sefidmazgi, A. G. 2020.  Causality Analysis of Climate and Ecosystem Time Series.  In _Evaluating Climate Change Impacts_ , 139–162. Chapman and Hall/CRC.
  * Shimizu et al. (2006) Shimizu, S.; Hoyer, P. O.; Hyvärinen, A.; Kerminen, A.; and Jordan, M. 2006.  A linear non-Gaussian acyclic model for causal discovery.  _Journal of Machine Learning Research_ , 7(10).
  * Siddiqi et al. (2022) Siddiqi, S. H.; Kording, K. P.; Parvizi, J.; and Fox, M. D. 2022.  Causal mapping of human brain function.  _Nature reviews neuroscience_ , 23(6): 361–375.
  * Smirnov (1939) Smirnov, N. V. 1939.  On the estimation of the discrepancy between empirical curves of distribution for two independent samples.  _Bull. Math. Univ. Moscou_ , 2(2): 3–14.
  * Smith et al. (2011) Smith, S. M.; Miller, K. L.; Salimi-Khorshidi, G.; Webster, M.; Beckmann, C. F.; Nichols, T. E.; Ramsey, J. D.; and Woolrich, M. W. 2011.  Network modelling methods for FMRI.  _Neuroimage_ , 54(2): 875–891.
  * Spirtes, Glymour, and Scheines (2000) Spirtes, P.; Glymour, C. N.; and Scheines, R. 2000.  _Causation, prediction, and search_.  MIT press.
  * Taylor (1998) Taylor, W. G. 1998.  _Effect of El Nino Southern Oscillation (ENSO) on British Columbia and Yukon winter weather_.  Aquatic and Atmospheric Science Division, Pacific and Yukon Region ….
  * Trifunov, Shadaydeh, and Denzler (2022) Trifunov, V.-T.; Shadaydeh, M.; and Denzler, J. 2022.  Time Series Causal Link Estimation under Hidden Confounding using Knockoff Interventions.  In _A causal view on dynamical systems, NeurIPS 2022 workshop_.
  * Wahl, Ninad, and Runge (2023) Wahl, J.; Ninad, U.; and Runge, J. 2023.  Vector causal inference between two groups of variables.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 37, 12305–12312.
  * Yao et al. (2020) Yao, L.; Chu, Z.; Li, S.; Li, Y.; Gao, J.; and Zhang, A. 2020.  A Survey on Causal Inference.  arXiv:2002.02770.
  * Zscheischler, Janzing, and Zhang (2012) Zscheischler, J.; Janzing, D.; and Zhang, K. 2012.  Testing whether linear equations are causal: A free probability theory approach.  _arXiv preprint arXiv:1202.3779_.