# Emerging Statistical Machine Learning Techniques for Extreme Temperature
Forecasting in U.S. Cities

Kameron B. Kinast
School of Mathematics and Statistics
Rochester Institute of Technology
Rochester, New York 14623
<EMAIL_ADDRESS>
&Ernest Fokoué
School of Mathematics and Statistics
Rochester Institute of Technology
Rochester, New York 14623
<EMAIL_ADDRESS>

###### Abstract

In this paper, we present a comprehensive analysis of extreme temperature
patterns using emerging statistical machine learning techniques. Our research
focuses on exploring and comparing the effectiveness of various statistical
models for climate time series forecasting. The models considered include
Auto-Regressive Integrated Moving Average, Exponential Smoothing, Multilayer
Perceptrons, and Gaussian Processes. We apply these methods to climate time
series data from five most populated U.S. cities, utilizing Python and Julia
to demonstrate the role of statistical computing in understanding climate
change and its impacts. Our findings highlight the differences between the
statistical methods and identify Multilayer Perceptrons as the most effective
approach. Additionally, we project extreme temperatures using this best-
performing method, up to 2030, and examine whether the temperature changes are
greater than zero, thereby testing a hypothesis.

_Keywords_ time series $\cdot$ weather forecasting $\cdot$ statistical models
$\cdot$ climate science

##  1 Introduction

For many decades, climate scientists have been concerned about the impact of
climate change. Climate change increases the risk to both natural and human
systems, and the degree of risk depends on various factors, including extreme
hot temperatures [1]. The 2015 Paris Agreement was established to limit the
increase in temperatures to 1.5∘C above pre-industrial levels (1850-1900) [2].
Achieving this goal would reduce exposure to climate-related risks, such as
heatwaves, droughts, and extreme precipitation events [1].

Understanding, modeling, and forecasting weather patterns pose ongoing
challenges. Time series forecasting is a technique used to predict future
observations by analyzing past values. Scientists employ various machine
learning and traditional approaches to analyze and predict these events. For
example, Kumar and Middey [3] used a hybrid of random forest and
autoregressive integrated moving average (ARIMA) model to project extreme
climate indicators. The objective of this paper is to provide a comprehensive
review of the effectiveness of commonly used classical and machine learning
methods for time series forecasting.

Numerous comparative studies have compared classical and machine learning
methods for time series forecasting. For instance, Hill and other colleagues
compared the neural network model with classical models such as the Box-
Jenkins model, single exponential smoothing model, and naive model, and
claimed that the neural network model performed the best[4]. Nesreen and other
authors [5] conducted a comparative study of different machine learning models
for time series forecasting and found that multilayer perceptron (MLP) and
Gaussian processes (GP) regression were the best models for forecasting with
M3 competition data or different types of time series data. Since classical
methods still hold significance in time series analysis and forecasting, this
paper will study four different time series forecasting methods:

  1. 1.

Auto-Regressive Integrated Moving Average

  2. 2.

Exponential Smoothing

  3. 3.

Multilayer Perceptron

  4. 4.

Gaussian Processes

As mentioned, the Paris Agreement aims to limit the temperature increase to
1.5 degrees Celsius to mitigate the vulnerability to severe climate effects.
This research paper investigates the rate of temperature change over a 28-year
period (2002-2030). Hypothetically, this paper performs time series
forecasting to test the hypothesis that the rate of temperature change is
greater than 0 ($\beta_{1}>0$). It is important to note that this research
paper analyzes a specific 28-year period, unlike the IPCC’s statement, which
expresses concerns about temperature increase over the pre-industrial period
spanning more than 150 years.

This paper explores various aspects of traditional approaches and machine
learning methods for time series forecasting using climate time series data.
The research will utilize Python [6] and Julia [7] to demonstrate the role of
statistical computing in understanding climate change. First, an overview of
the essential characteristics of the proposed statistical methods is provided.
Next, these methods are applied to climate time series data, and the
predictive performance of the proposed models is summarized. Finally, the
optimal statistical model is used to forecast temperatures for the next seven
years, up to 2030.

##  2 Statistical Models

###  2.1 Auto-Regressive Integrated Moving Average

The Auto-Regressive Integrated Moving Average (ARIMA) model is a combination
of three model components: the autoregressive (AR) model, the moving average
(MA) model, and the integrated method. An autoregressive model of order p,
denoted as AR(p), can be expressed as [8]:

$X_{t}=\phi_{1}X_{t-1}+\phi_{2}X_{t-2}+\ldots+\phi_{p}X_{t-p}+\epsilon_{t}$
(1)

Here, $X_{t}$ represents the value of the time series, which is defined as a
function of past values $X_{t-1},X_{t-2}\ldots X_{t-p}$, along with the
constant terms $\phi_{p}$ and the error term $\epsilon_{t}$. On the other
hand, a moving average model of order q, denoted as MA(q), can be expressed
as:

$X_{t}=\theta_{1}\epsilon_{t-1}+\theta_{2}\epsilon_{t-2}+\ldots+\theta_{q}\epsilon_{t-q}+\epsilon_{t}$
(2)

In this case, the time series is defined as a function of the errors of past
predictions $\epsilon_{t-1},\epsilon_{t-2},\ldots,\epsilon_{t-q}$, with
$\theta_{q}$ representing the constant terms. The integrated method is used to
deal with non-stationary time series through differencing. The ARIMA model is
typically denoted as ARIMA(p,d,q), where p is the order of the autoregressive
model, q is the order of the moving average model, and d is the order of
differencing in the integrated method.

The ARIMA model is based on the Box-Jenkins methodology, which involves an
iterative three-step process to develop a forecasting model for the time
series [9]:

  1. 1.

Ensure the time series is stationary. The Augmented Dickey-Fuller (ADF) test
can determine whether the time series is stationary or not.

  2. 2.

Identify the model using the autocorrelation function (ACF) and the partial
autocorrelation function (PACF). The ACF assesses the correlation between time
series considering all the lags or intervals between time periods, while the
PACF only considers specific lags.

  3. 3.

Estimate the model parameters using goodness-of-fit tests. Model selection is
usually based on the Akaike information criterion (AIC) and the Bayesian
information criterion (BIC), defined as follows:

$AIC=2k-2\ln(L)\mbox{ and }BIC=k\ln(n)-2\ln(L)$ (3)

The AIC measures the information value of the model using maximum likelihood
estimates, denoted as $L$, and the number of parameters in the model, denoted
as $k$. The BIC is similar to the AIC but includes a larger penalty term that
takes into account the number of observations, denoted as $n$, in the data.

###  2.2 Exponential Smoothing

Exponential Smoothing (ETS) is a traditional time series approach that
forecasts future values based on weighted averages of past observations,
assuming that current observations have more weight than past observations. By
applying the Holt-Winters additive method, the ETS model with the observed
time series $X_{t}$ can be expressed as [10]:

$\displaystyle
Level:L_{t}=\alpha(X_{t}-S_{t-p})+(1-\alpha)(L_{t-1}+T_{t-1})\hskip
28.45274pt0<\alpha<1$ (4) $\displaystyle
Trend:T_{t}=\beta(L_{t}-L_{t-1})+(1-\beta)T_{t-1}\hskip 28.45274pt0<\beta<1$
(5) $\displaystyle Seasonal:S_{t}=\gamma(X_{t}-L_{t})+(1-\gamma)S_{t-p}\hskip
28.45274pt0<\gamma<1$ (6)

Here, $p$ represents the number of periods in a seasonal cycle (e.g.,
quarterly, monthly, weekly), and $\alpha$, $\beta$, and $\gamma$ are smoothing
parameters. The forecasting equation at time t is given by:

$F_{t+n}=L_{t}+nT_{t}+S_{t-p+n},\hskip 28.45274ptn=1,2,3,....$ (7)

In this equation, $n$ represents the number of time periods ahead of the
current time period $t$.

###  2.3 Multilayer Perceptrons

The Multilayer Perceptron (MLP) is a widely used type of feedforward
artificial neural network in machine learning applications. The training
algorithm for an MLP is known as backpropagation, which involves adjusting the
network’s weights to minimize the difference between predicted outputs and the
actual outputs of the training data. The structure typically consists of an
input layer, one or more hidden layers of densely connected neurons with non-
linear activation functions, and an output layer. The input layer usually
consists of lagged observations of weather variables, while the output layer
produces forecasts of future weather conditions.

In this study, the MLP model’s structure consists of three layers: the input
layer, hidden layer, and output layer. The input layer has three neurons
representing three past observations in a sequence, and the output neuron
represents the forecasted observation. Let ${\bf x}\in\mathbb{R}^{p}$ be the
input vector, ${\bf z}^{(1)}\in\mathbb{R}^{m_{1}}$ be the vector of hidden,
latent variables for the first hidden layer, and $y\in\mathbb{R}$ be the
output. Assuming the MLP has three neurons in the input layer, four hidden
neurons in the hidden layer, and one neuron in the output layer, we have:

$\displaystyle{\bf z}^{(1)}$ $\displaystyle=$ $\displaystyle\sigma({\bf
W}^{(1)}{\bf x}+{\bf b}^{(1)})$ (8) $\displaystyle y$ $\displaystyle=$
$\displaystyle{\bf W}^{(2)}{\bf z}^{(1)}+{\bf b}^{(2)},$

Here, $\sigma(\cdot)$ is a nonlinear activation function applied element-wise
to the output of the linear transformations ${\bf W}^{(1)}{\bf x}+{\bf
b}^{(1)}$ and ${\bf W}^{(2)}{\bf z}^{(1)}+{\bf b}^{(2)}$, where ${\bf
W}^{(1)}\in\mathbb{R}^{m_{1}\times p}$ and ${\bf
b}^{(1)}\in\mathbb{R}^{m_{1}}$ are the weight matrix and bias vector for the
first hidden layer, and ${\bf W}^{(2)}\in\mathbb{R}^{m_{2}\times m_{1}}$ and
${\bf b}^{(2)}\in\mathbb{R}^{m_{2}}$ are the weight matrix and bias vector for
the output layer. In this case, the ReLU (Rectified Linear Unit) activation
function is used, which sets any negative values to 0 (Equation 9). This
function is computationally efficient and provides a simple nonlinear
transformation.

$\sigma(u)=max(0,u)$ (9)

To prevent overfitting in the MLP model, a regularization technique called
dropout is applied. Dropout randomly ignores selected neurons during training,
reducing the model’s sensitivity to specific neuron weights. The common
dropout rate used is $p=0.2$ in each hidden layer. With the dropout rate
represented as $d$, the hidden layer from Equation 8 can be described as [11]:

$\displaystyle d^{(1)}\sim\mbox{Bernoulli}(p),$ (10)
$\displaystyle\hat{z}^{(1)}=d^{(1)}\cdot z^{(1)}$ $\displaystyle
z^{(2)}=\phi^{(2)}(\Sigma W^{(2)}\hat{z}^{(1)}+b^{(2)}).$

This study utilizes the Adaptive Moment Estimation (Adam) optimization method
as the optimizer for the MLP model. To estimate the model’s losses, a loss
function is included, allowing the updating of parameters to minimize the
losses. Since the prediction involves real values, this MLP model is
considered a regression predictive modeling problem. Therefore, the default
loss function, Mean Squared Error (MSE), is used.

###  2.4 Gaussian Processes

Gaussian Processes (GPs) are powerful tools for time series analysis, offering
a flexible and non-parametric approach to model complex and non-linear
relationships in data. Fundamentally, a GP is a collection of random variables
in which all finite-dimensional distributions are joint Gaussian distributions
for any finite number [12]. In the context of time series analysis, a GP can
be viewed as a probability distribution over the set of all possible time
series functions, with the mean and covariance functions defining the
properties of the GP.

The goal of GP is to learn the underlying distribution from the training data,
and in order to do that, GP uses the method of Bayesian inference. Assuming
that the underlying function generating the time series values at time $t$ is
denoted by $f_{t}$, GP regression models the function $f_{t}$ as a Gaussian
process with a mean function $\mu_{t}$ and covariance function $\Sigma_{t}$.
Suppose the time series ${X_{t}}$ is given the observed values up to time $t$,
i.e., $\mathcal{D}_{t}={X_{1},X_{2},\ldots,X_{t}}$, the conditional
distribution of $X_{t+1}$ given $\mathcal{D}_{t}$ is then given by:

The goal of GP is to learn the underlying distribution from the training data,
achieved through the method of Bayesian inference. Assuming that the
underlying function generating the time series values at time $t$ is denoted
as $f_{t}$, GP regression models the function $f_{t}$ as a Gaussian process
with a mean function $\mu_{t}$ and a covariance function $\Sigma_{t}$. Given
the observed values of the time series ${X_{t}}$ up to time $t$, denoted as
$\mathcal{D}_{t}={X_{1},X_{2},\ldots,X_{t}}$, the conditional distribution of
$X_{t+1}$ given $\mathcal{D}_{t}$ is given by:

$\displaystyle p(X_{t+1}\mid\mathcal{D}_{t})$ $\displaystyle=$
$\displaystyle\mathcal{N}(\mu_{t+1},\Sigma_{t+1}),$ (11)

The posterior mean and covariance are computed as:

$\displaystyle\mu_{t+1}$ $\displaystyle=$ $\displaystyle
k_{t+1}^{\top}(K_{t}+\sigma^{2}I)^{-1}X_{t},$ (12) $\displaystyle\Sigma_{t+1}$
$\displaystyle=$ $\displaystyle
k(X_{t+1},X_{t+1})-k_{t+1}^{\top}(K_{t}+\sigma^{2}I)^{-1}k_{t+1},$ (13)

Here,
$k_{t+1}=[k(X_{1},X_{t+1}),k(X_{2},X_{t+1}),\ldots,k(X_{t},X_{t+1})]^{\top}$,
$K_{t}$ is the $t\times t$ matrix with entries $K_{i,j}=k(X_{i},X_{j})$,
$\sigma^{2}$ is the noise variance, and $I$ is the identity matrix. The kernel
functions $k(X_{i},X_{j})$ specify the correlation between the predicted and
observed time series values.

GPs often employ a kernel composition approach where multiple kernels are
combined through addition or multiplication to shape the resulting
distribution. In this study, the GP kernel composition is expressed as:

$k(X_{i},X_{j})=k_{Noise}(X_{i},X_{j})+k_{Periodic}(X_{i},X_{j})+k_{RBF}(X_{i},X_{j})$
(14)

where $k_{Noise}(X_{i},X_{j})$ represents the noise kernel capturing the noise
aspects of the time series, $k_{Periodic}(X_{i},X_{j})$ represents the
periodic kernel capturing seasonal patterns, and $k_{RBF}(X_{i},X_{j})$
represents the RBF (Radial Basis Function) kernel introducing non-linear
trends.

The $k_{Periodic}(X_{i},X_{j})$ and $k_{RBF}(X_{i},X_{j})$ are expressed as:

$\displaystyle k_{Periodic}(X_{i},X_{j})$ $\displaystyle=$
$\displaystyle\exp\left(-\frac{2\sin^{2}\left(\pi|X_{i}-X_{j}|/p\right)}{\ell^{2}}\right)$
(15)

and

$\displaystyle k_{RBF}(X_{i},X_{j})$ $\displaystyle=$
$\displaystyle\exp\left(-\frac{(X_{i}-X_{j})^{2}}{2\ell^{2}}\right).$ (16)

###  2.5 Model Selection Criterion

To determine the optimal time series forecasting model, the performance is
evaluated using testing datasets. The performance of the four different time
series forecasting models is assessed using the following statistical
measures: root mean squared error (RMSE) (Equation 17) and mean absolute error
(MAE) (Equation 18).

$RMSE=\sqrt{\frac{\sum_{i=1}^{N}(y_{i}-\hat{y}_{i})^{2}}{N}}$ (17)

$MAE=\frac{1}{N}\sum_{i=1}^{N}|y-\hat{y}|$ (18)

In these equations, $y_{i}$ represents the actual observed values,
$\hat{y}_{i}$ represents the predicted values, and $n$ represents the number
of observations in the testing dataset. These statistical measures provide a
quantitative assessment of the accuracy and performance of the forecasting
models by capturing the differences between the predicted and actual values.

##  3 Application to Climate Time Series Data

###  3.1 Data used and study area

The climate data utilized in this study is obtained from the National Centers
for Environmental Information’s (NCEI’s) Climate Data Online (CDO), a
repository managed by the National Oceanographic and Atmospheric
Administration (NOAA). CDO provides comprehensive summaries of historical
daily land surface observations from locations worldwide.

For this analysis, the climate dataset focuses on two key time series
measurements: maximum temperatures and minimum temperatures. Maximum and
minimum temperatures are recorded in Fahrenheit. The dataset consists of
observations collected from January 1st, 1950, to December 31st, 2022,
encompassing data from five airport sites across the United States. The
selected airport sites, representing various regions of the United States, are
as follows:

  1. 1.

Houston

  2. 2.

Chicago

  3. 3.

Boston

  4. 4.

San Francisco

  5. 5.

Miami

Given the extensive size of the climate dataset, which spans nearly 75 years
of data points for each variable, a mini-batch approach is employed to focus
on the most recent twenty years of data. Consequently, the dataset is divided
into a training set and a testing set in an 80% to 20% ratio, allowing for a
comparative analysis of the time series forecasting methods. With this split,
the training set comprises approximately 16 years of time series observations,
while the testing set encompasses the most recent four years, specifically
from 2019 to 2022.

###  3.2 ARIMA Results

The ACF and PACF are critical statistical measures for analyzing time series
data. For each dataset, plots of the ACF and PACF are examined to determine
the appropriate model order. The selection of the best models is based on the
Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC),
and the results are summarized in Table 1.

Variables | City | AIC score | BIC score | Model
---|---|---|---|---
TMAX | Houston | 3.38 | 20558.26 | ARIMA(88,0,0)
Chicago | 4.04 | 24473.82 | ARIMA(109,0,0)
Boston | 4.06 | 24590.41 | ARIMA(108,0,0)
San Francisco | 3.08 | 18412.60 | ARIMA(19,0,0)
Miami | 2.30 | 14700.06 | ARIMA(100,0,0)
TMIN | Houston | 3.28 | 20070.57 | ARIMA(101,0,0)
Chicago | 3.66 | 22267.02 | ARIMA(106,0,0)
Boston | 3.23 | 19856.57 | ARIMA(108,0,0)
San Francisco | 1.95 | 12845.68 | ARIMA(97,0,0)
Miami | 2.55 | 16098.39 | ARIMA(106,0,0)
Table 1: Model Order Selection for ARIMA model of Maximum and Minimum
Temperatures among five cities.

###  3.3 ETS Results

The analysis of the time series data revealed the presence of seasonal
components in the maximum and minimum temperatures. Since the training set
spans 16 years, it suggests the presence of a recurring pattern within each
year. Accordingly, the ETS model, also known as Triple Exponential Smoothing
(TES), is employed without considering a trend but with a periodicity of 16
seasons.

###  3.4 MLP Results

Prior to training the MLP model, the training set is split into sequences of
three observations, which serve as input neurons for the model. Both the
training and testing datasets are divided into batches of 12 samples to update
the model’s parameters iteratively. The model undergoes 2,000 epochs of
training to optimize its parameters and minimize the loss function.

###  3.5 GP Results

For the GP model in this study, the kernel composition includes the following
components:

  * •

a noise kernel with a standard deviation of 1.0

  * •

a periodic kernel with a length scale of 0.0, a standard deviation of 1.0, and
a periodicity of 1.0

  * •

and an RBF kernel with a length scale of 4.0 and standard deviation of 0.0.

The parameter values for each kernel were determined through multiple trial-
and-error iterations using the root mean squared error (RMSE) as the
criterion, as summarized in Table 2. With this fixed kernel composition, the
GP model utilizes a zero mean function since it is the standard mean function
for time series observation data, implying that the mean value is consistently
zero.

Kernel Composition | RMSE Values
---|---
N(1) + P(0,1,1) + RBF(4,0) | 12.51
N(1) + P(0,1,1) | 12.81
P(0,1,1) + RBF(4,0) | 12.81
N(1) + P(0,1,2) + RBF(4,0) | Failed due to infinite sin values
N(1) + P(0,1,0) + RBF(4,0) | Failed due to infinite sin values
N(1) + P(0,1,1) + RBF(2,0) | 12.80
N(1) + P(0,1,1) + RBF(5,0) | 12.83
N(1) + P(0,1,1) + RBF(4,4) | 76
N(1) + P(0,1,0) + RBF(4,4) | 76
N(1) + P(0,1,1) + RBF(3,0) | 12.86
N(1) + P(1,1,1) + RBF(4,0) | Failed due to infinite sin values
Table 2: Trials and Errors with Values for Kernel in Gaussian Process, using
Houston’s Maximum Temperatures

###  3.6 Comparative Results

The forecasting values of the four methods, along with the testing set of
observed values from the NCEI’s CDO, are presented in Figure 1 for maximum
temperatures and Figure 2 for minimum temperatures. The graphs illustrate that
the predicted values by the MLP model closely match the actual values, unlike
the predictions from the other three models.

Table 3 and Table 4 display the model performance of the four methods. Upon
examination of the tables, it is observed that the fixed kernel composition of
the GP model fails to produce valid predictions for the maximum and minimum
temperatures of three cities due to infinite values in the sine function. GP
models are known for their instability and ill-conditioning, which means that
the fixed kernel composition might be ineffective for these cities, leading to
errors in inverse matrices.

Based on the tables, the RMSE and MAE measures indicate that the MLP model
outperforms the other three methods. Hence, it can be concluded that the MLP
model demonstrates the best performance in forecasting the maximum and minimum
temperatures for these five cities. The overall ranking of the four models,
based on their performance, is as follows: MLP, ARIMA, ETS, and GP. However,
if the model performance is assessed solely based on these two cities, the GP
model performs better than the ETS model.

Model | Measures | Houston | Chicago | Boston | San Francisco | Miami
---|---|---|---|---|---|---
ARIMA | RMSE | 11.24 | 16.75 | 13.74 | 8.62 | 5.45
MAE | 8.90 | 14.10 | 11.33 | 6.69 | 4.29
ETS | RMSE | 28.14 | 30.28 | 25.54 | 14.62 | 7.11
MAE | 25.63 | 25.09 | 20.87 | 12.07 | 5.92
MLP | RMSE | 6.09 | 8.24 | 8.25 | 4.99 | 3.54
MAE | 4.42 | 6.33 | 6.53 | 3.81 | 2.29
GP | RMSE | 12.51 | – | 17.94 | – | –
MAE | 10.36 | – | 15.24 | – | –
Table 3: Model Performance Measures for Forecasting on Maximum Temperatures among five cities. Model | Measures | Houston | Chicago | Boston | San Francisco | Miami
---|---|---|---|---|---|---
ARIMA | RMSE | 11.60 | 16.75 | 11.59 | 5.50 | 6.82
MAE | 9.65 | 12.90 | 9.43 | 4.53 | 5.48
ETS | RMSE | 21.50 | 23.67 | 22.52 | 9.41 | 7.32
MAE | 18.65 | 19.29 | 18.48 | 8.00 | 5.56
MLP | RMSE | 5.71 | 6.56 | 5.42 | 2.81 | 4.00
MAE | 4.02 | 5.01 | 4.07 | 2.13 | 2.80
GP | RMSE | 13.77 | – | 16.36 | – | –
MAE | 11.96 | – | 13.91 | – | –
Table 4: Model Performance Measures for Forecasting on Minimum Temperatures
among five cities.

(a) Houston

(b) Chicago

(c) Boston

(d) San Francisco

(e) Miami

Figure 1: Testing set of maximum temperatures values and fitted values
predicted by the four methods.

(a) Houston

(b) Chicago

(c) Boston

(d) San Francisco

(e) Miami

Figure 2: Testing set of minimum temperatures values and fitted values
predicted by the four methods.

###  3.7 Temperature Projections

Based on the analysis of the figures and the evaluation of model performance,
it is evident that the MLP model outperforms the other three statistical
models. Figure 3 presents the projections of maximum temperatures, while
Figure 4 illustrates the projections of minimum temperatures for all five
cities from 2023 to 2030, utilizing the MLP model. The observations values
from 2002 to 2022 are incorporated in the model.

Each graph in Figures 3 and 4 includes a regression equation displayed in the
bottom left corner. Notably, all regression equations demonstrate that the
coefficient $\beta$ is greater than 0, indicating a positive rate of change in
the temperatures.

(a) Houston

(b) Chicago

(c) Boston

(d) San Francisco

(e) Miami

Figure 3: Observation values of the maximum temperatures of five cities with
the forecast values using MLP method from 2003 to 2030.

(a) Houston

(b) Chicago

(c) Boston

(d) San Francisco

(e) Miami

Figure 4: Observation values of the minimum temperatures of five cities with
the forecast values using MLP method from 2003 to 2030.

##  4 Conclusion

The foundation of the developed time series methods is based on recent
advancements in machine learning and statistical techniques, which provide an
optimal framework for each weather forecasting model. In this paper, we
explored the feasibility of four different statistical models for time series
forecasting of temperatures in five cities. Each city has a unique temperature
range, influenced by factors such as its location (coastal or mountainous). We
observed that temperatures in all cities exhibit significant seasonal
components but no trend.

Based on statistical measures such as RMSE and MAE, the MLP model demonstrates
the highest accuracy in predicting maximum and minimum temperatures across all
cities. The ARIMA model performs second best, followed by ETS and GP. While
this paper suggests that the MLP model is the most effective for time series
forecasting of temperatures, it is important to consider certain limitations
in model development, such as the focus on only one structure of the MLP model
and a fixed kernel composition of GP. Additionally, due to the large dataset
size, running these methods on the entire dataset can cause crashes on typical
computers. Therefore, this research focuses on a mini-batch of a 20-year
dataset instead of 75 years. With access to high-performance computing
technology, these machine learning techniques could potentially deliver even
better model performance.

The MLP model, as the best model selection, forecasts the temperatures of all
five cities up to 2030, and all forecasting performances confirm that the
hypothesis is true where there is a positive rate of change, $\beta>0$.
Therefore, it indicates that the temperatures, either minimum or maximum, are
increasing within only 28 years (2002-2030). However, as presented in the
results, the MLP model is not able to forecast extreme temperatures.

It is worth noting that time series analysis encompasses numerous methods
beyond the scope of this research, and statistical techniques continue to
advance. Furthermore, there are areas for further research expansion. For
instance, an empirical comparison study could be extended to include
additional machine learning methods for weather time series forecasting, such
as Long-Short Term Memory (LSTM), Support Vector Machine (SVM), and
Convolutional Neural Network (CNN).

As part of this research, I have created a GitHub repository where all of the
code related to my research is stored. The repository is publicly accessible
and intended to serve as a resource for anyone interested in the work I have
done. By making my code available to the public, I hope to encourage others to
build upon my research and to contribute to the wider community of scholars in
this field. You can find the repository at
https://github.com/kamkinast24/Time-Series-Approach.

## References

  * Masson-Delmotte et al. [2018] V. Masson-Delmotte, P. Zhai, H.-O. Pörtner, D. Roberts, J. Skea, P.R. Shukla, A. Pirani, W. Moufouma-Okia, C. Péan, S. Connors R. Pidcock, J.B.R. Matthews, Y. Chen, X. Zhou, M.I. Gomis, E. Lonnoy, T. Maycock, M. Tignor, and T. Waterfield (eds.).  Global warming of 1.5°c. an ipcc special report on the impacts of global warming of 1.5°c above pre-industrial levels and related global greenhouse gas emission pathways, in the context of strengthening the global response to the threat of climate change, sustainable development, and efforts to eradicate poverty.  _Summary for Policymakers_ , pages 3–24, 2018.  URL https://doi.org/10.1017/9781009157940.001.
  * Rogelj et al. [2016] Joeri Rogelj, Michel den Elzen, Niklas Höhne, Taryn Fransen, Hanna Fekete, Harald Winkler, Roberto Schaeffer, Fu Sha, Keywan Riahi, and Malte Meinshausen.  Paris agreement climate proposals need a boost to keep warming well below 2°c.  _Nature_ , 534(7609):631–639, 2016.  doi:10.1038/nature18307.  URL https://doi.org/10.1038/nature18307.
  * Kumar and Middey [2023] N. Kumar and A. Middey.  Extreme climate index estimation and projection in association with enviro-meteorological parameters using random forest-arima hybrid model over the vidarbha region, india.  _Environmental Monitoring and Assessment_ , 195(380), 2023\.  doi:10.1007/s10661-022-10902-2.  URL https://doi.org/10.1007/s10661-022-10902-2.
  * Hill et al. [1996] Timothy Hill, Marcus O’Connor, and William Remus.  Neural network models for time series forecasts.  _Management Science_ , 42:1082–1092, 11 1996.  doi:10.1287/mnsc.42.7.1082.  URL https://doi.org/10.1287/mnsc.42.7.1082.
  * Ahmed et al. [2010] Nesreen K. Ahmed, Amir F. Atiya, Neamat El Gayar, and Hisham El-Shishiny.  An empirical comparison of machine learning models for time series forecasting.  _Econometric Reviews_ , 29(5-6):594–621, 2010\.  doi:10.1080/07474938.2010.481556.  URL https://doi.org/10.1080/07474938.2010.481556.
  * Van Rossum and Drake Jr [1995] Guido Van Rossum and Fred L Drake Jr.  _Python reference manual_.  Centrum voor Wiskunde en Informatica Amsterdam, 1995.
  * Bezanson et al. [2012] Jeff Bezanson, Stefan Karpinski, Viral B Shah, and Alan Edelman.  Julia: A fast dynamic language for technical computing.  _arXiv preprint arXiv:1209.5145_ , 2012.
  * Shumway and Stoffer [2017] Robert H. Shumway and David S. Stoffer.  _ARIMA Models_ , pages 75–163.  Springer International Publishing, Cham, 2017.  ISBN 978-3-319-52452-8.  doi:10.1007/978-3-319-52452-8_3.  URL https://doi.org/10.1007/978-3-319-52452-8_3.
  * Dimri et al. [2020] Tripti Dimri, Shamshad Ahmad, and Mohammad Sharif.  Time series analysis of climate variables using seasonal arima approach.  _Journal of Earth System Science_ , 129(1), 06 2020.  URL https://doi.org/10.1007/s12040-020-01408-x.
  * Srivastava et al. [2016] Prashant K. Srivastava, Tanvir Islam, Sudhir K. Singh, George P. Petropoulos, Manika Gupta, and Qiang Dai.  Forecasting arabian sea level rise using exponential smoothing state space models and arima from topex and jason satellite radar altimeter data.  _Meteorological Applications_ , 23(4):633–639, Dec 2016.  doi:10.1002/met.1585.  URL https://doi.org/10.1002/met.1585.
  * Srivastava et al. [2014] Nitish Srivastava, Geoffrey Hinton, Alex Krizhevsky, Ilya Sutskever, and Ruslan Salakhutdinov.  Dropout: A simple way to prevent neural networks from overfitting.  _Journal of Machine Learning Research_ , 15(56):1929–1958, 2014.  URL http://jmlr.org/papers/v15/srivastava14a.html.
  * Quadrianto et al. [2010] Novi Quadrianto, Kristian Kersting, and Zhao Xu.  _Gaussian Process_ , pages 428–439.  Springer US, Boston, MA, 2010.  ISBN 978-0-387-30164-8.  doi:10.1007/978-0-387-30164-8_324.  URL https://doi.org/10.1007/978-0-387-30164-8_324.
  * Lewinson [2020] Eryk Lewinson.  _Multilayer perceptrons for time series forecasting_ , chapter 10, page 361–376.  Packt Publishing Ltd., 2020.


11institutetext: University of Science - Faculty of Mathematics and Computer
Science
Ho Chi Minh City, Vietnam 22institutetext: Vietnam National University Ho Chi
Minh City, Vietnam 33institutetext: SimulaMet, Oslo, Norway

# Correlation visualization under missing values:
a comparison between imputation and direct parameter estimation methods

Nhat-Hao Pham1,2,∗ Khanh-Linh Vo1,2,∗ Mai Anh Vu1,2 Thu Nguyen3
Michael A. Riegler3 Pål Halvorsen3 and Binh T. Nguyen1,2
1University of Science - Faculty of Mathematics and Computer Science
Ho Chi Minh City  Vietnam
2 Vietnam National University Ho Chi Minh City  Vietnam
3SimulaMet - Dept. of Holistic Systems  Oslo  Norway
111denotes equal contribution Nhat-Hao Pham⋆ 1 122 Khanh-Linh Vo⋆ 1122 Mai Anh
Vu denotes equal contribution 1 122 Thu Nguyen 33 Michael A. Riegler 33 Pål
Halvorsen 33 Binh T. Nguyen 1122

###### Abstract

Correlation matrix visualization is essential for understanding the
relationships between variables in a dataset, but missing data can seriously
affect this important data visualization tool. In this paper, we compare the
effects of various missing data methods on the correlation plot, focusing on
two randomly missing data and monotone missing data. We aim to provide
practical strategies and recommendations for researchers and practitioners in
creating and analyzing the correlation plot under missing data. Our
experimental results suggest that while imputation is commonly used for
missing data, using imputed data for plotting the correlation matrix may lead
to a significantly misleading inference of the relation between the features.
In addition, the most accurate technique for computing a correlation matrix
(in terms of RMSE) does not always give the correlation plots that most
resemble the one based on complete data (the ground truth). We recommend using
DPER [1], a direct parameter estimation approach, for plotting the correlation
matrix based on its performance in the experiments.

###### Keywords:

missing data correlation plot data visualization.

##  1 Introduction

A correlation plot is a powerful and important visualization tool to summarise
datasets, reveal relationships between variables, and identify patterns in the
given dataset. However, missing data is a common situation in practice. There
are two main approaches to computing correlation when the data contains
missing entries. The first is to impute missing data and then estimate the
correlation from imputed data. The second way is to estimate the correlation
matrix directly from missing data.

Various imputation methods have been proposed to handle missing data. For
example, some traditional methods that have been used for missing data include
mean imputation and mode imputation. In addition, multivariate imputation by
chained equations (MICE) [2] is an iterative approach that can provide the
uncertainty of the estimates. Furthermore, machine learning techniques train
on the available data to learn patterns and relationships and then are used to
predict the missing values. There are many typical approaches, including tree-
based models [3], matrix completion [4], or clustering-based methods such as
KNN Imputation (KNNI). Finally, deep learning is an emerging trend for
handling missing data, due to the excellent prediction of deep learning
architectures. Notable deep learning methods for imputation can be Generative
Adversarial Imputation Nets (GAIN) [5] and Graph Imputer Neural Network (GINN)
[6]. These methods have achieved high accuracy but require significant
computational resources and are best suited for data-rich cases.

For the direct parameter estimate approach, a simple way to estimate the
covariance matrix from the data is to use case deletion. However, the accuracy
of such an approach may be low. Recently, direct parameter estimation
approaches were proposed to estimate the parameters from the data directly.
For example, the EPEM algorithm [7] estimates these parameters from the data
by maximum likelihood estimation for monotone missing data. Next, the DPER
algorithm [1] is based on maximum likelihood estimation for randomly missing
data [1]. The computation for the terms in the covariance matrix is based on
each pair of features. Next, the PMF algorithm [8] for missing dataset with
many observed features by using principle component analysis on the portion of
observed features and then using maximum likelihood estimation to estimate the
parameters. Such algorithms are scalable and have promising implementation
speed.

Since the correlation plot is an important tool in data analysis, it is
necessary to understand the impact of missing data on the correlation plot
[9]. This motivates us to study the effects of various imputation and
parameter estimation techniques on the correlation plot. In summary, our
contributions are as follows. First, we examine the effect of different
imputation and direct parameter estimation methods on the correlation matrix
and the correlation matrix heat maps for randomly missing data and monotone
missing data. Second, we show that imputation techniques that yield the lowest
RMSE may not produce the correlation plots that most resemble the one based on
completed data (the ground truth). Third, we illustrate that relying solely on
RMSE to choose an imputation method and then plotting the correlation heat map
based on the imputed data from that imputation method can be misleading. Last
but not least, our analysis process using the Local RMSE differences heatmaps
and Local difference (matrix subtraction) for correlation provides
practitioners a strategy to draw inferences from correlation plots under
missing data.

##  2 Methods Under Comparison

This section details techniques under comparison. These include some state-of-
the-art techniques (DPER, GAIN, GINN, SoftImpute, ImputePCA, missForest), and
some traditional but remaining widely used methods (KNNI, MICE, Mean
imputation). More details are as follows, (Except for DPER [1], which directly
estimates the covariance matrix from the data, the remaining methods involve
imputing missing values prior to calculating the covariance matrix.)

  * •

Mean imputation imputes the missing entries for each feature by the mean of
the observed entries in the corresponding feature.

  * •

Multiple Imputation by Chained Equations (MICE) [10] works by using a series
of regression models that fit the observed data. In each iteration of the
process, the missing values for a specific feature are imputed based on the
other features in the dataset.

  * •

MissForest [3] is an iterative imputation method based on a random forest. It
averages over many unpruned regression or classification trees, allowing a
multiple imputation scheme and allowing estimation of imputation error.

  * •

K-nearest neighbor imputation (KNNI) is based on the K-nearest neighbor
algorithm. For each sample to be imputed, KNNI finds a set of K-nearest
neighbors and imputes the missing values with a mean of the neighbors for
continuous data or the mode for categorical data.

  * •

SoftImpute [4] is an efficient algorithm for large matrix factorization and
completion. It uses convex relaxation to produce a sequence of regularized
low-rank solutions for completing a large matrix.

  * •

ImputePCA (iterative PCA) [11] minimizes the criterion, and the imputation of
missing values is achieved during the estimation process. Since this algorithm
is based on the PCA model, it takes into account the similarities between
individuals and the relationships between variables.

  * •

Expectation Maximization algorithm (EM) [12] is a maximum likelihood approach
from incomplete data. It iteratively imputes missing values and updates the
model parameters until convergence. The algorithm consists of two steps:
E-step estimates missing data given observed data, and M-step estimates model
parameters given complete data.

  * •

Generative Adversarial Imputation Networks (GAIN) [5] adapts the GAN
architecture where the generator imputes the missing data, while the
discriminator distinguishes between observed and imputed data. The generator
maximizes the discriminator’s error rate, while the discriminator minimizes
classification loss.

  * •

Graph Imputation Neural Networks (GINN) [6] is an autoencoder architecture
that completes graphs with missing values. It consists of a graph encoder and
a graph decoder, which utilize the graph structure to propagate information
between connected nodes for effective imputation.

  * •

Direct Parameter Estimation for Randomly missing data (DPER algorithm) [1] is
a recently introduced algorithm that directly finds the maximum likelihood
estimates for an arbitrary one-class/multiple-class randomly missing dataset
based on some mild assumptions of the data. Since the correlation plot does
not consider the label, the DPER algorithm used for the experiments is the
DPER algorithm for one class dataset.

##  3 Experiments

###  3.1 Experiment Settings

This section details the experimental setup for comparing the effects of the
methods presented in Section 2 on correlation plots with missing data. For the
experiments on randomly missing data, we used the Iris and Digits datasets and
generated randomly missing data with rates from $10\%$ to $50\%$. Here, the
missing rate is the ratio between the number of missing entries and the total
number of entries. For the monotone missing pattern, we removed a $50\%$ piece
of the image in the right corner with $40\%,50\%$, and $60\%$ height and width
in the MNIST dataset [13]. For MNIST, as KNNI and GINN [6] cannot produce the
results within 3 hours of running, we exclude them from the plot of results.
For each dataset, we normalize the data and run the experiments on a Linux
machine with 16 GB of RAM and 4 cores, Intel i5-7200U, 3.100GHz.

Dataset | # Classes | # Features | # Samples
---|---|---|---
Iris | 3 | 4 | 150
Digits | 10 | 64 | 1797
MNIST [13] | 10 | 784 | 70000
Table 1: The description of datasets used in the experiments.

To throughout investigate the effects of missing data on the correlation
heatmap, we employ three types of heatmap plots:

  * •

Correlation Heatmaps illustrates the correlation matrix and employ a blue-to-
white-to-red color gradient to illustrate correlation coefficients. The
colormap utilizes distinct colors: blue, white, and red, which are represented
by the values $[(0,0,1)$, $(1,1,1)$, $(1,0,0)]$ respectively, as defined by
the color map of the matplotlib library. Blue indicates negative values, red
denotes positive values, and white indicates median (0) values. The colormap’s
range is bounded within vmin = -1 and vmax = 1, which are parameters of
matplotlib that linearly map the colors in the colormap from data values vmin
to vmax.

  * •

Local RMSE Difference Heatmaps for Correlation portray the RMSE differences on
a cell-by-cell basis, comparing correlation matrices produced by each
estimation technique with the ground truth correlation matrix. We adopt a
different color gradient, which is a green-to-white color map with values
$[(1,1,1)$, $(0,0.5,0)]$, further enhancing the clarity of these
visualizations. The intensity of the green colors corresponds to the magnitude
of the difference between the correlation matrix and the ground truth
(depicted as white).

  * •

Local Difference (Matrix Subtraction) Heatmaps for Correlation show the local
differences for a cell-by-cell comparison (matrix subtraction) between
correlation matrices generated from each estimation technique and the ground
truth correlation matrix. This plot shows whether the difference between the
two plots is positive or negative, providing a clear indication of the
direction of the difference.

For each dataset, we provide five heatmap plots. The first two plots display
the method outcomes across varying levels of missing value, with the first
showing the Correlation Heatmaps and the second showing the Local RMSE
Difference Heatmaps. The two subsequent plots provide a more detailed
examination, specifically highlighting the highest missing rate of the
Correlation Heatmaps and Local RMSE difference heatmaps. This enables a more
thorough analysis of color variations that might have been challenging to
detect in the first two plots due to spatial limitations. The final plot, also
focusing on the highest missing rate, shows the Local Difference (Matrix
Subtraction) Heatmaps for Correlation, clearly indicating the direction of the
difference.

We use gray to differentiate null values, especially when all cells in a
feature share the same values, often near the image edges. This gray shade
distinguishes null values from 0 values, which are depicted in white.
Furthermore, the plot arrangement follows a descending order of RMSE, focusing
on the highest missing rate’s RMSE due to its significant impact on color
intensity across rates. To maintain consistency with the correlation plot, we
represent the ground truth heatmaps difference in white. In addition, in the
Local RMSE Difference Heatmaps for Correlation charts, we integrate the
calculation of the RMSE, which represents the difference between the two
correlation matrices while excluding null positions.

Moreover, we include the dense rank in ascending order alongside it. Dense
ranking is a method for giving ranks to a group of numbers. The ranks are
given so that there are no gaps between them. The next rank isn’t skipped if
some numbers share the same rank. For instance, let’s take this set of values:
[89, 72, 72, 65, 94, 89, 72]. The resulting dense ranks would be: [1, 2, 2, 2,
3, 3, 4]. Below each Local RMSE Difference Heatmaps plot, we denote the rank
and the RMSE value using the following structure:
$(<\text{rank}>)\text{RMSE:}<\text{RMSE value}>$.

###  3.2 Result & Analysis

Figure 1: The RMSE difference of the correlation matrix between each method
and the ground truth is depicted across missing rates ranging from $10\%\text{
to }50\%$ for (a) Iris and Digits dataset with randomly missing patterns and
(b) MNIST dataset with monotone missing patterns.

####  3.2.1 Randomly missing data

(a) Correlation Heatmaps

(b) Local RMSE Difference Heatmaps for Correlation

Figure 2: Heatmaps for the Iris dataset across missing rates from $10\%$ to
$50\%$

A clear trend emerges in the RMSE values across the Iris and Digits datasets
as depicted in the initial two line plots in Figure 1. As the missing rate
increases, there is a corresponding rise in the RMSE difference. This trend is
also observed in the Local RMSE Difference Heatmaps for Correlation at 3(b)
5(b), where the color intensity of the green shades increases in both Figure
2(b) and Figure 4(b). The charts also show that the RMSE rankings remain
consistent by comparing the changing ranks within the Local RMSE difference
heatmaps in the last plot of Figures 2 and 4. This indicates that the findings
about different correlation matrix estimation methods stay the same across
different missing rates.

For the Iris Dataset, at its highest missing rate in the Local RMSE Difference
Heatmaps for Correlation plot in Figure 3(b), with Mean Imputation as the
baseline, it’s worth noting that EM and MissForest also exhibit RMSE values
that are close to the baseline. Seven of these methods perform better than the
baseline, while MissForest and GINN demonstrate higher RMSE values (0.2429 and
0.3416, respectively). Notably, DPER, MICE, and GAIN stand out, as they
exhibit RMSE values lower than the median (RMSE = 0.1383, represented by
KNNI). These methods are highly recommended for estimating the correlation
matrix in this dataset.

The visualization of the Local Difference (Matrix Subtraction) Heatmaps for
Correlation in Figure 3(c) indicates that the discrepancies are both positive
and negative, represented by a mixed color between blue and white. However, we
observe that the cells with an intense green color for GINN, MissForest, Mean
Impute, and EM in Figure 3(b) appear blue in Figure 3(c). This implies that
features with greater divergence compared to the ground truth tend to have a
negative tendency in their differences, suggesting an underestimation of the
correlation between these features.

(a) Correlation Heatmaps

(b) Local RMSE Difference Heatmaps for Correlation

(c) Local Difference (Matrix Subtraction) Heatmaps for Correlation

Figure 3: Heatmaps for the Iris dataset at a missing rate of $50\%$

For the Digits dataset, in the first row of the plot layout, several methods
outperform Mean Impute and have RMSE better than the median value of 0.097,
namely KNN, MICE, SoftImpute, and DPER. Interestingly, while EM has a better
RMSE than Mean Impute, their RMSE values are quite similar. Like the Iris
dataset, EM, Mean Impute, and MissForest exhibit nearly identical RMSE values,
while GINN still displays higher RMSE values. Both ImputePCA and GAIN show
larger RMSE values compared to the baseline.

(a) Correlation Heatmaps

(b) Local RMSE Difference Heatmaps for Correlation

Figure 4: Heatmaps for the Digits dataset across missing rates from $10\%$ to
$50\%$

(a) Correlation Heatmaps

(b) Local RMSE Difference Heatmaps

(c) Local Difference (Matrix Subtraction) Heatmaps for Correlation

Figure 5: Heatmaps for the Digits dataset at a missing rate of $50\%$

Visualizing the Local Difference (Matrix Subtraction) Heatmaps for Correlation
in Figure 5(c) reveals both positive and negative disparities. An intriguing
pattern emerges: certain lines running parallel to the diagonal consistently
display intense green colors in Figure 5(b), suggesting significant
differences. These areas of intense green indicate that these features tend to
underestimate correlation, as evidenced by the corresponding positions of
these cells in Figure 5(c) appearing as blue. These positions correspond to
cells with very high correlation values, as visible in Figure 5(a). This
phenomenon is consistent across SoftImpute, EM, Mean Impute, MissForest, GINN,
and GAIN.

Regarding correlation estimation, methods like KNNI, MICE, SoftImpute, and
DPER perform better than the median RMSE value and are thus recommended.
However, practitioners should exercise caution when using SoftImpute, as it
may underestimate the correlation of some highly correlated features.

Throughout the analysis above, techniques including KNN, MICE, SoftImpute, and
DPER consistently demonstrate their effectiveness in estimating correlation
matrices for the Digits dataset. Within the context of the Iris dataset, DPER,
MICE, and GAIN emerge as reliable methods. Moreover, our visual observations
highlight certain methods’ inclination to underestimate highly correlated
features in the correlation matrix within the Digits dataset.

####  3.2.2 Monotone missing data

In contrast to the Iris and Digits datasets, the results for the plots in
Figure 6(b) across varying missing rates show minimal variation among methods.
Most methods, except for GAIN, exhibit closely aligned RMSE values, which
demonstrate the highest RMSE and are visually apparent in the RMSE plot. All
of these methods, except for GAIN, are recommended for estimating the
correlation matrix for this dataset.

(a) Correlation Heatmaps

(b) Local RMSE Difference Heatmaps for Correlation

Figure 6: Heatmaps for the MNIST dataset across missing rates from $40\%$ to
$50\%$

(a) Correlation Heatmaps

(b) Local RMSE Difference Heatmaps for Correlation

(c) Local Difference (Matrix Subtraction) Heatmaps for Correlation

Figure 7: Heatmaps for the MNIST dataset at a missing rate of $50\%$

####  3.2.3 Analysis

In summary, from the observations for all the datasets, we see that DPER and
MICE consistently provide plots that most resemble the ground truth.

From the experiments, one can see that a low RMSE may not imply a correlation
plot reflecting the true relation between the features. For example, at a 50%
missing rate on the Digits dataset, even though softImpute and Missforest have
lower RMSE than DPER, DPER still clearly resembles the ground truth plot. To
understand why, suppose that we have a dataset of two features where the
covariance matrix is

$\mathbf{\Sigma}=\begin{pmatrix}\sigma_{1}&\sigma_{12}\\\
\sigma_{12}&\sigma_{2}\end{pmatrix}.$

Then
$RMSE=\sqrt{(\hat{\sigma}_{1}-\sigma_{1})^{2}+(\hat{\sigma}_{2}-\sigma_{2})^{2}+2(\hat{\sigma}_{12}-\sigma_{12})^{2}},$
and whether
$((\hat{\sigma}_{1}-\sigma_{1})^{2},(\hat{\sigma}_{2}-\sigma_{2})^{2},(\hat{\sigma}_{12}-\sigma_{12})^{2})$
receives the values $(0.1,0.5,0.1)$ or $(0.2,0.2,0.2)$ will result in the same
values of RMSE.

Also, it is important to emphasize that more than relying solely on RMSE
values for selecting an imputation method, and then using the imputed data
based on that method to create correlation plots and draw inferences can lead
to misleading conclusions. Considering Local RMSE differences in heatmaps and
Local difference (matrix subtraction), heatmaps are vital to ensure precise
method selection.

##  4 Conclusions

In this work, we examined the effects of missing data on correlation plots.
Crucially, we found that the techniques yielding the lowest RMSE do not
consistently produce correlation plots that closely resemble those based on
complete data (the ground truth). This could be because the RMSE is an average
measure. Meanwhile, human eyes capture local differences in the matrix, while
visual perception is better at capturing local differences, which could be
better visualized by representing the local square difference on each cell of
the heatmap (employed by the Local RMSE Difference Heatmaps for Correlation in
our experiments). In addition, the ground truth data itself may contain noise.
Therefore, relying solely on RMSE values for selecting an imputation method
and then using the imputed data based on that method to create a correlation
plot and draw inferences can lead to misleading conclusions. Considering Local
RMSE differences in heatmaps and Local difference (matrix subtraction),
heatmaps are vital to ensure proper inferences from correlation plots when
there is missing data.

Moreover, our analysis consistently highlights the efficiency of correlation
estimation methods such as MICE and DPER across all three datasets. Notably,
DPER is scalable and consistently produces the correlation plot that best
resembles the ground truth for both small and large datasets. We hypothesize
that such a performance is related to the asymptotic consistency property of
DPER and will investigate this theoretically in the future.

## References

  * [1] Thu Nguyen, Khoi Minh Nguyen-Duy, Duy Ho Minh Nguyen, Binh T. Nguyen, and Bruce Alan Wade.  Dper: Direct parameter estimation for randomly missing data.  Knowledge-Based Systems, 240:108082, 2022.
  * [2] S van Buuren and Karin Groothuis-Oudshoorn.  mice: Multivariate imputation by chained equations in r.  Journal of statistical software, pages 1–68, 2010.
  * [3] Daniel J Stekhoven and Peter Bühlmann.  Missforest-non-parametric missing value imputation for mixed-type data.  Bioinformatics, 28(1):112–118, 2012.
  * [4] Rahul Mazumder, Trevor Hastie, and Robert Tibshirani.  Spectral regularization algorithms for learning large incomplete matrices.  The Journal of Machine Learning Research, 11:2287–2322, 2010.
  * [5] Jinsung Yoon, James Jordon, and Mihaela van der Schaar.  GAIN: missing data imputation using generative adversarial nets.  CoRR, abs/1806.02920, 2018.
  * [6] Indro Spinelli, Simone Scardapane, and Aurelio Uncini.  Missing data imputation with adversarially-trained graph convolutional networks.  Neural Networks, 129:249–260, sep 2020.
  * [7] Thu Nguyen, Duy H.M. Nguyen, Huy Nguyen, Binh T. Nguyen, and Bruce A. Wade.  Epem: Efficient parameter estimation for multiple class monotone missing data.  Information Sciences, 567:1–22, 2021.
  * [8] Thu Nguyen, Nhan Thanh Phan, Ha Van Hoang, Pål Halvorsen, Michael A Riegler, and Binh T Nguyen.  Pmf: Efficient parameter estimation for data sets with missing data in some features.  Available at SSRN 4260235.
  * [9] Matthias Kraus, Katrin Angerbauer, Juri Buchmüller, Daniel Schweitzer, Daniel A Keim, Michael Sedlmair, and Johannes Fuchs.  Assessing 2d and 3d heatmaps for comparative analysis: An empirical study.  In Proceedings of the 2020 CHI Conference on Human Factors in Computing Systems, pages 1–14, 2020.
  * [10] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay.  Scikit-learn: Machine learning in Python.  Journal of Machine Learning Research, 12:2825–2830, 2011.
  * [11] Julie Josse and François Husson.  missmda: a package for handling missing values in multivariate data analysis.  Journal of statistical software, 70:1–31, 2016.
  * [12] Arthur P Dempster, Nan M Laird, and Donald B Rubin.  Maximum likelihood from incomplete data via the em algorithm.  Journal of the Royal Statistical Society: Series B (Methodological), 39(1):1–38, 1977.
  * [13] Yann LeCun.  The mnist database of handwritten digits.  http://yann. lecun. com/exdb/mnist/, 1998.