2021

[2]Juliane Müller

1]Artificial Intelligenc Team, Teikametrics, Boston, MA 02210, United States

2]Computational Research Division, Lawrence Berkeley National Laboratory, 1
Cyclotron Rd, Berkeley, CA 94720, United States

3]Earth and Environmental Sciences Area, Lawrence Berkeley National
Laboratory, Calvin Rd, Berkeley, CA 94705, United States

4]Water Security Research Group, International Institute for Applied Systems
Analysis, Schloßpl. 1, 2361 Laxenburg, Austria

# Long-Term Missing Value Imputation for Time Series Data Using Deep Neural
Networks

Jangho Park111This work was done when Jangho Park was with Lawrence Berkeley
National Laboratory<EMAIL_ADDRESS>Bhavna Arora  Boris Faybishenko
Gilberto Pastorello  Charuleka Varadharajan  Reetik Sahu  Deborah Agarwal  [
[  [  [

###### Abstract

We present an approach that uses a deep learning model, in particular, a
MultiLayer Perceptron (MLP), for estimating the missing values of a variable
in multivariate time series data. We focus on filling a long continuous gap
(e.g., multiple months of missing daily observations) rather than on
individual randomly missing observations. Our proposed gap filling algorithm
uses an automated method for determining the optimal MLP model architecture,
thus allowing for optimal prediction performance for the given time series. We
tested our approach by filling gaps of various lengths (three months to three
years) in three environmental datasets with different time series
characteristics, namely daily groundwater levels, daily soil moisture, and
hourly Net Ecosystem Exchange. We compared the accuracy of the gap-filled
values obtained with our approach to the widely-used R-based time series gap
filling methods ImputeTS and mtsdi. The results indicate that using an MLP for
filling a large gap leads to better results, especially when the data behave
nonlinearly. Thus, our approach enables the use of datasets that have a large
gap in one variable, which is common in many long-term environmental
monitoring observations.

###### keywords:

Missing value imputation, Environmental data, Machine learning, Hyperparameter
optimization, Derivative-free optimization, Surrogate models

##  1 Introduction and Literature Review

Time series data are recorded in diverse application areas ranging from earth
sciences to healthcare, finance, traffic, etc. Time series often have missing
values (gaps) due to, for example, sensor failures, collection errors, or lack
of resources pedro2009 ; yozgatligil2013comparison ; kalteh2009imputation ;
aissia2017multivariate . However, often complete time series datasets are
required for analysis or when these data are used as inputs to numerical
models. For environmental science applications, such time series gaps could
coincide with crucial times where hydrological or biogeochemical fluxes or
rates are high. As an example, high stream stage conditions or fluctuations in
temperature can promote carbon and nutrient cycling in hyporheic zone
environments gu2012riparian ; arora2019evaluating . Thus, it becomes crucial
to impute missing time series values with highly accurate estimates. There are
two categories of time series missing value imputation methods, namely
multivariate and univariate methods thi2019edtwbi ; moritz2015comparison . The
first class estimates missing values by exploiting relationships between
variables. The second class solely relies on available observations of the
time series whose gaps must be filled. In this paper, we study a special case
where a large number of consecutive values are missing in a single variable of
a multivariate time series.

Many different methods have been used to deal with different characteristics
of missing values. For example, for small gaps in univariate time series
(individual missing values), the last observation before the gap can be
carried forward or the next observation after the gap backward to estimate
missing values. Other approaches, which are implemented in ImputeTS
moritz2017imputets  provide linear, spline interpolation, or Kalman smoothing
methods. Machine learning-based algorithms have also been studied for gap
filling purposes. For example, 9237768  estimated missing values in univariate
CO2 concentrations and air temperature data. This approach, however, is not
able to take advantage of multidimensional data.

For the imputation of missing values in the multivariate case, the $K$-nearest
neighbor method has been successfully applied for imputing data in the field
of medical science batista2002study . However, this and other such
applications have been limited to imputation problems where values are missing
only at random. Recently, recurrent neural networks che2018recurrent ;
NEURIPS2018_734e6bfc  and generative adversarial networks luo2018multivariate
; zhang2021missing  have been used for imputing values in health-care and
public air quality datasets. These studies focused exclusively on the
performance and accuracy of the proposed approaches for filling individual
randomly missing values.

MLPs have also been successfully used for imputing missing time series values
across areas of differing complexity and applications. For example,
lingras2008evolutionary ; zhong2007rationalizing  showed that time delayed
deep neural network models can impute missing values in univariate hourly
traffic volumes. The referenced works focus on filling relatively short gaps
(up to 168 data points) and a genetic algorithm was employed for optimizing
the hyperparameters of the used gap filling methods. mahabbati2020comparison
used an MLP to fill missing values in 16 meteorological drivers and three key
ecosystem fluxes in the regional Australian and New Zealand flux tower
network. The authors analyzed the performance of gap filling methods by
manually changing the lags (the number of historical time steps used for
making the next time step prediction). To estimate the missing flux values,
they used multiple meteorological drivers such as air temperature, wind speed,
soil water content, ground heat flux, and net radiation. The architecture of
their MLP was determined by using a grid search method provided in scikit-
learn sklearn . In another study, kim2020gap  focused on filling the gaps in
methane flux measurements using multiple input variables including net
ecosystem exchange (NEE) of CO2, latent heat flux, sensible heat flux, global
radiation, outgoing long wave radiation, air temperature, soil temperature,
relative humidity, vapor pressure deficit, friction velocity, air pressure,
and water table height. They used an MLP with two hidden layers and determined
the number of nodes by using a three-fold cross‐validation method.

Unlike the referenced studies, the main contribution of our work is a method
for filling a long continuous gaps (e.g., multiple continuous months of
missing daily observations) in a single variable of a multivariate time series
dataset. Another key difference between the previous studies and our research
is the size of the datasets used. Our focus is on smaller datasets that have
fewer variables than the cited works. Usually, these types of datasets are
hard to use and tend to be discarded in analysis and modeling efforts. Our
approach would enable greater use of time series datasets that have multi-year
gaps in a variable, which are common in many long-term environmental
monitoring observations. In contrast to prior approaches using MLPs, our
proposed gap filling method optimizes the MLP model architecture and the
length of the lag with an efficient method that is based on surrogate models.
Our method requires fewer resources because it is automated and uses adaptive
sampling rather than evolutionary or trial-and-error approaches.

In this paper, we present a gap filling approach that is based on the concept
of Time-Delayed Deep Neural Networks (TDNN) waibel1989phoneme  with the goal
of filling a long continuous gap in one variable of a multivariate time
series. The concept of TDNN is suitable for handling temporal time series that
are a characteristic feature of environmental datasets. This research extends
our previous work muller2019surrogate  in which we used deep neural networks
to predict groundwater levels for multiple years into the future using
supporting variables such as temperature, precipitation, and river flux data.
In the cited work, we optimized the network architectures with an automated
hyperparameter tuning method, which we adopt and modify in the present work to
enable gap filling. We assume that the target time series variable has one
large gap and that the supporting variables that explain the target variable
are fully observed. Our previous work showed that MLPs are successful in
making accurate predictions, and therefore, we use an MLP for gap filling in
the present study. Our research will have an impact on application areas where
gap-free long-term observations are particularly useful.

The remainder of this paper is organized as follows. Section 2 provides a
brief description of MLP models and hyperparameter optimization (HPO) used for
gap filling. Section 3 presents the results of our numerical study in which we
compare the performance of our proposed method with ImputeTS
moritz2017imputets  and mtsdi junger2018package  for filling gaps of various
lengths (three months to three years) of three different environmental
datasets. Finally, Section 4 concludes the paper and provides future research
directions.

##  2 Multilayer Perceptrons and Hyperparameter Optimization for Gap Filling

In this section, we provide brief descriptions of MLPs and the hyperparameter
optimization method used to determine its architecture.

###  2.1 Multilayer Perceptron (MLP)

In this paper, we use MLPs for regression purposes to map input to output
features. An MLP is a feed-forward artificial neural network. A generic MLP
consists of an input layer, at least one hidden layer, and an output layer
with each layer containing multiple nodes. Each node in the input and output
layers corresponds to an input feature and target variable, respectively. In
order to construct an input layer that encodes temporal information, we
structure each input such that it contains multiple time steps. We define the
amount of historical information in an input as lag $l$. Thus, with $T$ time
periods, an input is defined as $\\{x_{t-l},\ldots,x_{t-1},x_{t}\\}$, where
$t=l+1,\ldots,T$. Our goal is to predict the value at the next time step, and
thus the output layer has only a single node, which corresponds to the target
variable $Y$ that needs to be gap-filled. With this architecture, the MLP is
trained by minimizing over a set of weights and biases using a predefined loss
function such as a root mean squared error or other metric that reflects how
well the MLP approximates the data. The root mean squared error is the square
root of the average difference between the actual (ground truth) and predicted
values. In the first step of this training, linear combinations of input nodes
with initial weights are constructed. Each linear combination is passed to all
nodes in the first hidden layer. An activation function in the nodes receives
this value and returns a transformed value. The transformed values from all
nodes in the first hidden layer are given to the nodes in the second hidden
layer for further transformation. This process repeats until the output layer
returns the final transformed value. The loss function to be minimized
reflects the error between the transformed output values and the true values.
By iteratively adjusting the weights and biases associated with each node, the
loss function is minimized. We use the Rectified Linear Unit (ReLU) activation
function, $f(\cdot)=\max(0,\cdot)$ nair2010rectified  and our training
objective function (loss) is the mean squared error (MSE). The iterative
optimization of the weights and biases is performed with the Adam optimizer
kingma2014adam , which is a type of backpropagation algorithm based on
gradient descent and the chain rule rumelhart1986learning .

Other essential hyperparameters that influence the performance of the MLP,
such as the batch size (size of batched samples in the gradient update), the
number of epochs (number of iterations used in training), the dropout rate
srivastava2014dropout , the number of hidden layers and nodes in each layer,
and the lag value $l$, are determined by the automated optimization method
described in the next section. We assume that the learning and decay rates of
the stochastic gradient descent method used for training the MLP are known and
fixed. For further details about MLPs, we refer the reader to
goodfellow2016deep .

###  2.2 Optimization of MLP Architectures for Gap Filling

In order to improve the predictive performance of MLPs (and deep learning
models in general), the hyperparameters (model architecture) should be
carefully chosen feurer2019hyperparameter . Different hyperparameter
optimization (HPO) methods have been developed and used in the literature,
including random search bergstra2012random , a combination of grid search and
trial and error larochelle2007empirical , Bayesian optimization
snoek2012practical , and genetic algorithms xie2017genetic .

In this work, we adapt the HPO method proposed in muller2019surrogate , which
interprets the training of a given architecture as an expensive black-box
function evaluation of the architecture’s performance. Thus, the goal is to
find the best hyperparameters within as few trials as possible by employing
computationally efficient methods. To this end, surrogate models (in
particular radial basis functions and Gaussian process models) are used to map
the hyperparameters to their respective performance. An initial experimental
design with $n_{0}$ different hyperparameter sets is created for which the
performance is evaluated by training the respective models. This initial
design is used to create the computationally cheap surrogate models. In each
iteration of the optimization, an auxiliary optimization problem is solved on
the computationally fast to evaluate surrogate model to determine the next
hyperparameters to be tried. The surrogate model is updated each time a new
hyperparameter set has been evaluated.

In order to adapt the HPO algorithm described in muller2019surrogate  for our
gap filling purposes, we only have to change the way we split the training and
validation datasets by taking into account the fact that some data are
missing. Since we include the lag (the number of historical time steps for
predicting the value at the next time step) in our list of hyperparameters to
be optimized, we have to create a lagged table of the data for each set of
hyperparameters we want to evaluate.

We denote the observation matrix with $N$ variables and time length $S$ as
$X=(x_{1},x_{2},\ldots,x_{S})^{T}\in\mathbb{R}^{S\times N}$, where
$x_{s}=(x_{s}^{1},x_{s}^{2},\ldots,x_{s}^{N})\in\mathbb{R}^{N}$ is an
observation vector at time $s,1\leq s\leq S$. Each $x_{s}^{d}$ is an
observation of the $d^{th}$ variable, where $d=1,\ldots,N$. We assume that the
first column in the matrix $X$ corresponds to the target vector $Y$ that
contains missing values, and the remaining columns are the supporting
variables. In the following, we denote $x^{d}_{s}=\text{NaN}$ if the value is
missing. For example, assume we have nine observations ($S=9$), one target
variable, and two supporting variables ($N=3$), then the observation matrix
and target value vector are, respectively,

$X=\blockarray{cccc}&123\\\ \block{c(ccc)}1x^{1}_{1}x^{2}_{1}x^{3}_{1}\\\
2x^{1}_{2}x^{2}_{2}x^{3}_{2}\\\ 3x^{1}_{3}x^{2}_{3}x^{3}_{3}\\\
4x^{1}_{4}x^{2}_{4}x^{3}_{4}\\\ 5\text{NaN}x^{2}_{5}x^{3}_{5}\\\
6\text{NaN}x^{2}_{6}x^{3}_{6}\\\ 7x^{1}_{7}x^{2}_{7}x^{3}_{7}\\\
8x^{1}_{8}x^{2}_{8}x^{3}_{8}\\\ 9x^{1}_{9}x^{2}_{9}x^{3}_{9}\\\ \hskip
5.69046pt,\hskip 28.45274ptY=\blockarray{cc}\block{c(c)}1&x^{1}_{1}\\\
2x^{1}_{2}\\\ 3x^{1}_{3}\\\ 4x^{1}_{4}\\\ 5\text{NaN}\\\ 6\text{NaN}\\\
7x^{1}_{7}\\\ 8x^{1}_{8}\\\ 9x^{1}_{9}\\\ \ .$

Here, the values for $x^{1}_{5}$ and $x^{1}_{6}$ are missing and the goal is
to estimate their values using an optimized MLP. Then, the lagged table for,
e.g., lag $l=1$ becomes

$\text{Input}=\blockarray{ccccccc}&123456\\\
\block{c(cccccc)}1x^{1}_{1}x^{2}_{1}x^{3}_{1}x^{1}_{2}x^{2}_{2}x^{3}_{2}\\\
2x^{1}_{2}x^{2}_{2}x^{3}_{2}x^{1}_{3}x^{2}_{3}x^{3}_{3}\\\
3x^{1}_{3}x^{2}_{3}x^{3}_{3}x^{1}_{4}x^{2}_{4}x^{3}_{4}\\\
4x^{1}_{4}x^{2}_{4}x^{3}_{4}\text{NaN}x^{2}_{5}x^{3}_{5}\\\
5\text{NaN}x^{2}_{5}x^{3}_{5}\text{NaN}x^{2}_{6}x^{3}_{6}\\\
6\text{NaN}x^{2}_{6}x^{3}_{6}x^{1}_{7}x^{2}_{7}x^{3}_{7}\\\
7x^{1}_{7}x^{2}_{7}x^{3}_{7}x^{1}_{8}x^{2}_{8}x^{3}_{8}\\\ \hskip
5.69046pt,\hskip
28.45274pt\text{Output}=\blockarray{cc}\block{c(c)}1&x^{1}_{3}\\\
2x^{1}_{4}\\\ 3\text{NaN}\\\ 4\text{NaN}\\\ 5x^{1}_{7}\\\ 6x^{1}_{8}\\\
7x^{1}_{9}\\\ .$

After the lagged tables have been constructed, we determine if they contain
any missing values (NaN). Any row that contains at least one NaN entry in the
input or output can be removed because the rows in the lagged table are
independent. Thus, for the above example, we obtain the following reduced
inputs and outputs

$\text{Reduced Input}=\blockarray{ccccccc}&123456\\\
\block{c(cccccc)}1x^{1}_{1}x^{2}_{1}x^{3}_{1}x^{1}_{2}x^{2}_{2}x^{3}_{2}\\\
2x^{1}_{2}x^{2}_{2}x^{3}_{2}x^{1}_{3}x^{2}_{3}x^{3}_{3}\\\
7x^{1}_{7}x^{2}_{7}x^{3}_{7}x^{1}_{8}x^{2}_{8}x^{3}_{8}\\\ ,\text{Reduced
Output}=\blockarray{cc}\block{c(c)}1&x^{1}_{3}\\\ 2x^{1}_{4}\\\ 7x^{1}_{9}\\\
\ .$

Once the reduced lagged input and output are obtained, they are divided into
training and validation datasets. For any architecture, the reduced lagged
training dataset is used to train the corresponding neural network model.
Since this training is done by stochastic gradient descent (and thus the final
trained model can be interpreted as the realization of a stochastic process),
we train the model $k$ times, each with a different random number seed. Each
of the $k$ trained models is validated against the reduced lagged validation
dataset, and the architecture’s performance is calculated as the sample
average over all $k$ values.

Finally, we impute the missing values sequentially starting from the first
missing value to the last missing value, using the previous estimated missing
value to impute the next missing value. In the small example above, we impute
the first missing value $x^{1}_{5}$ and denote it as $\hat{x}^{1}_{5}$ by
using $(x^{1}_{3},x^{2}_{3},x^{3}_{3},x^{1}_{4},x^{2}_{4},x^{3}_{4})$. Then,
the next missing value ($x_{6}^{1}$) is estimated based on $\hat{x}^{1}_{5}$,
i.e., $(x^{1}_{4},x^{2}_{4},x^{3}_{4},\hat{x}^{1}_{5},x^{2}_{5},x^{3}_{5})$.
In this way, all missing values are estimated.

##  3 Numerical Experiments

In this section, we demonstrate that our proposed method for gap filling of
multivariate time series data is general enough to be used in different
applications that have varying observation frequencies and data
characteristics. We describe three time series datasets for which we assume
that the target variable has one large gap that must be filled. These time
series include observed groundwater levels, simulated soil moisture, and
derived net ecosystem exchange from flux tower measurements. Below, we provide
details of the setup of the numerical experiments and alternative gap filling
algorithms that we use for comparison.

###  3.1 Description of the Datasets

####  3.1.1 Measured groundwater levels

Our first dataset contained groundwater level time series measured in a
groundwater monitoring well in Butte County in California (CA), United States.
The monitoring well 22222N01E28J001M, California Natural Resources Agency,
Periodic Groundwater Level Measurements,
https://data.cnra.ca.gov/dataset/periodic-groundwater-level-measurements is
200 m deep from the surface, and provides the groundwater elevation above the
mean sea level at daily time resolution. We obtained the data from the
California Natural Resources Agency. Our goal was to gap fill the time series
of groundwater levels using measurements of the supporting variables -
temperature, precipitation and river discharge. The nearest discharge
monitoring station (Butte Creek Durham) measures the daily discharge rate at
the Butte Creek about 8 km from the well333California Data Exchange Center,
Butte Creek Durham station,
https://cdec.water.ca.gov/webgis/?appid=cdecstation&sta=BCD. Temperature and
precipitation data were obtained from the Chico weather station located 7 km
from the well site444California Data Exchange Center, Chico Station,
https://cdec.water.ca.gov/webgis/?appid=cdecstation&sta=CHI. The dataset used
in the study consisted of approximately eight years (2010-2018) of daily
observations. For a more detailed description of the dataset and the data
sources, we refer the reader to muller2019surrogate ; sahu2020impact . Figure
1 shows the time series plots for all variables.

Figure 1: Time series plots for the groundwater level imputation test case.
Shown are (top to bottom) the daily groundwater levels (GW), temperature
(Temp), precipitation (Rain), and log-transformed riverflux.

####  3.1.2 Simulated soil moisture

Our second dataset was created with the HYDRUS-1D model hydrusmanual ;
vsimuunek2008modeling , a soil hydrologic simulation model. The HYDRUS-1D
model is a modular, freely-available, state-of-the-art, and widely-used soil
hydrologic model with many advanced and coupled water and reactive chemical
transport features including equivalent continuum and dual permeability
modeling approaches for preferential flow and transport. In particular, the
model simulates water, heat and solute transport in variably-saturated and
saturated porous media. The model has been extensively applied for modeling
across scales – from laboratory cores to watersheds arora2015integrated ;
baek2020assessment .

Here, we used the HYDRUS-1D model to simulate a representative soil column in
the Butte County in California to determine the changes in soil moisture
profiles over approximately seven years (2011-2018). As a result, our target
variable was the daily total water content (TWC) in the upper 80 cm of the
simulated soil column. The net amount of water entering the soil column from
the top (vTop) and the groundwater recharge were used as the supporting
variables for calculations. Figure 2 shows the time series plots for all
variables. Because simulation models usually do not lead to missing values, we
use it here as an ideal case study to investigate the general applicability of
our MLP-based gap filling approach.

Figure 2: Time series plots for the soil moisture imputation test case. Shown
are (top to bottom) the daily total water content (TWC), the net amount of
water entering the soil column from the top (vTop), and the groundwater
recharge.

####  3.1.3 Net ecosystem exchange data

Our third example application was the most recently produced FLUXNET dataset,
that is FLUXNET2015 dataset. It includes data on CO2, water, energy exchange,
and other meteorological and biological measurements pastorello2020fluxnet2015
; FLUXNET2015 . The eddy-covariance method is used to allow for the non-
destructive estimation of fluxes between atmosphere and biosphere for a single
site to global scale. FLUXNET datasets have been used in a wide range of
research areas ranging from soil microbiology to validation of large-scale
earth system models.

The target and supporting variables for the Morgan Monroe State Forest site,
Indiana, United States for 16 years (1999-2014) novick2016fluxnet2015  are
described in Table 1 and the time series are illustrated in Figure 3. Note
that none of the variables in the studied date range have missing values.

Table 1: Target and supporting variables with description of FLUXNET2015
dataset

| Variable | Units | Description
---|---|---|---
Target variable | NEE_VUT_USTAR50 |  $\mu$molCO2/m2s |  Net Ecosystem Exchange, using Variable Ustar Threshold (VUT) for each year, from 50 percentile of USTAR threshold
Supporting variables | TA_F |  ∘C |  Air temperature
| VPD_F | hecto Pascals |  Vapor Pressure Deficit
| SW_IN_F | Watt/m2 |  Shortwave radiation, incoming

Figure 3: Time series plots for the net ecosystem exchange test case. Shown
are (top to bottom) the hourly Net Ecosystem Exchange (NEE_VUT_USTAR50), Air
Temperature (TA_F), Vapor Pressure Deficit (VPD_F), and incoming shortwave
radiation (SW_IN_F).

In all three case studies, the supporting variables had been identified with
the help of domain science experts, but filter methods, recursive feature
elimination, or sensitivity analyses sahu2020impact  may prove useful to
further downselect the most important features. For all test cases, temporal
information such as the time of day, the day, and the month of the year the
measurements were taken were also included as additional supporting variables.
Also, all values of each feature were rescaled with the min-max normalization,
i.e., for each variable $x^{i}$, to transform all variables to have the same
scale. We found the minimum value $x^{i}_{\text{min}}$ and the maximum value
$x^{i}_{\text{max}}$, and we rescaled the values according to
$\tilde{x}^{i}=\frac{x^{i}-x^{i}_{\text{min}}}{x^{i}_{\text{max}}-x^{i}_{\text{min}}}$.

###  3.2 Setup of Numerical Experiments

####  3.2.1 HPO details

We implemented the HPO algorithm for gap filling described in Section 2.2 in
python ($version~{}3.7$) using PyTorch ($version~{}1.4.0$) paszke2017automatic
. All experiments were run on Ubuntu 16.04 with Intel® Xeon(R) CPU E3-1245 v6
@ 3.70GHz $\times$ 8, and 31.2 GB memory. We used the same parameter settings
as in muller2019surrogate  and we optimized over the following hyperparameters
and their ranges:

  * [noitemsep]

  * •

batch size $\in\\{10,15,\ldots,195,200\\}$,

  * •

epochs $\in\\{50,100,\ldots,450,500\\}$,

  * •

layers $\in\\{1,2,\ldots,5,6\\}$,

  * •

number of hidden nodes $\in\\{5,10,\ldots,45,50\\}$,

  * •

dropout rate $\in\\{0.0,0.1,\ldots,0.4,0.5\\}$,

  * •

lags $\in\\{30,35,\ldots,360,365\\}$.

The optimization hyperparameters and their possible values were the same for
all three test cases. Note that all parameters were mapped to consecutive
integers for the optimization, and there is a finite number of possible
combinations (architectures). With the provided options, there are 7,588,800
possible MLP architectures, and it is computationally intractable to train
that many networks, thus motivating the use of automated HPO.

We assumed that each hidden network layer has the same number of hidden nodes
and the same dropout rate. We use the MSE as performance measure in tuning the
hyperparematers and we split the data into 85% training data and 15%
validation data.

In the HPO algorithm, we set $n_{0}=10$ as the initial experimental design
size (we used 10 different hyperparameter sets and we trained the models for
each to obtain the model performance, which we use to initialize the surrogate
model). We stopped the algorithm after $n=50$ hyperparameter sets had been
tried. In order to take into account the stochasticity that arises from using
the stochastic gradient descent method for training the models, we computed
the hyperparameter performance as the sample average over multiple training
trials. For example, for the daily groundwater level application and the
HYDRUS-1D application, we trained the model for each hyperparameter set five
times using the same training dataset, and the performance was defined as the
average MSE over those trials. Due to the computational expense associated
with training the MLP models for the FLUXNET data (because of the large amount
of training data), we computed the performance only over one training trial
for each hyperparameter set.

####  3.2.2 Description of gaps in time series data

In order to test our developed gap filling method, we created an artificial
continuous gap in the target time series. We studied different instances in
which the gaps have varying lengths and are located at different times of the
year. By using artificial gaps, we have the ability to compare the gap-filled
values to the true held-out observations as well as analyze and compare the
performance of different gap filling methods. In all three test cases, the
supporting variables did not have gaps.

For the groundwater application, we created gaps of three, six and twelve
months of daily observations. For each case, we used 12 different problem
instances in which we removed the data of different months of the year, which
allowed us to analyze a possible dependence of the performance of the gap
filling methods on the trend characteristics of the missing data. For example,
for the three month gap, we assumed that groundwater levels were missing (i)
from January through March 2012, (ii) from February through April 2012, (iii)
from March through May 2012, etc. For the HYDRUS-1D application, we assumed
that the total water content values were missing for 12 and 15 months,
respectively. For the 12-month gap, we assumed that the total water content
values were missing (i) from January 2013 through December 2013, (ii) from
February 2013 through January 2014, and (iii) from March 2013 through February
2014. For the 15 months gap, we considered 9 cases where the missing values
were staring from the months of April, May, till December. Finally, for the
hourly FLUXNET data, we considered between 2,159 hours (3 months) and 35,063
hours (36 months) of missing values in the Net Ecosystem Exchange data
(NEE_VUT_USTAR50).

####  3.2.3 Gap filling algorithms used in the numerical comparison

We compared our proposed gap filling method with state-of-the-art imputation
packages, including ImputeTS moritz2017imputets  and mtsdi junger2018package .
ImputeTS is a univariate time series missing value imputation package
implemented in R. We used the na.seadec function with an interpolation
algorithm. The na.seadec function imputes missing values with deseasonalized
time series data and adds the seasonality again. It has the option to identify
the number of observations before the seasonal pattern repeats. However, in
our applications, this feature did not yield satisfactory results, and thus we
manually set the daily time series frequency, i.e., the frequency of 365.
mtsdi is an expectation-maximization algorithm based multivariate time series
missing value imputation package in R. We used the default settings, i.e., the
mnimput function with the spline method, and we set the spline smooth control
to 7. With mtsdi, we included indicator variables such as the month, week, and
day.

###  3.3 Results and Discussion

####  3.3.1 Groundwater level predictions

Table 2 presents a comparison of the root mean squared error (RMSE) computed
between the predicted values and the true (held-out) values for the
groundwater use cases with three, six, and twelve months of missing daily data
values, respectively. Shown are the RMSEs of our proposed MLP-based gap
filling method, ImputeTS, and mtsdi. The RMSE is given in meters (m). We show
the optimal hyperparameters obtained with our method when using the Gaussian
process (GP) and the radial basis function (RBF) as surrogate model,
respectively.

For three months of missing daily observations (top section of Table 2), our
results show that the lowest RMSE was achieved by different methods depending
on the date range for which the data are missing. For example, for the three
months that span November 2012-January 2013 (see also Figure 4), the
groundwater levels show a linear growth. All gap filling methods captured this
increasing trend and ImputeTS yielded the best performance in terms of RMSE.
This is due to the fact that ImputeTS uses na.seadec with an interpolation
algorithm. In other words, ImputeTS performed seasonally decomposed missing
value imputation by linear interpolation which yielded excellent agreement
with the held-out data for the missing data segment that had linear behavior.
Also the MLP-based gap filling methods were able to capture the linear trend,
but the predictions were not as accurate. mtsdi introduced strong oscillations
around the overall growing trend. Moreover, mtsdi performed overall the worst
for all date ranges for the three month test case.

Table 2: RMSE values (m) for the different gap filling methods for three, six,
and twelve months of missing values in daily groundwater levels and optimal
hyperparameters ([batch size, # epochs, # layers, # nodes per layer, dropout
rate, # lags]) for both RBF- and GP-based HPO.66footnotemark: 6

Missing date range (yyyy-mm-dd) | #obs | RMSE |  | Hyperparameters
---|---|---|---|---
From | To |  | MLP RBF | MLP GP | ImputeTS | mtsdi |  | MLP RBF | MLP GP
2012-01-01 | 2012-03-31 | 91 | 0.34 | 0.38 | 0.56 | 1.23 |  | [145, 150, 1, 45, 0, 330] | [190, 450, 4, 50, 0, 270]
2012-02-01 | 2012-04-30 | 90 | 0.45 | 0.56 | 0.49 | 1.18 |  | [120, 400, 1, 25, 0, 250] | [95, 450, 3, 40, 0, 210]
2012-03-01 | 2012-05-31 | 92 | 0.54 | 0.61 | 0.97 | 1.34 |  | [195, 400, 1, 45, 0, 250] | [170, 350, 2, 20, 0, 290]
2012-04-01 | 2012-06-30 | 91 | 0.56 | 0.62 | 0.72 | 1.38 |  | [175, 450, 1, 40, 0, 240] | [70, 500, 1, 25, 0, 345]
2012-05-01 | 2012-07-31 | 92 | 0.61 | 0.62 | 0.61 | 1.17 |  | [130, 450, 4, 45, 0, 330] | [170, 350, 1, 30, 0, 240]
2012-06-01 | 2012-08-31 | 92 | 0.42 | 0.41 | 0.39 | 1.64 |  | [170, 250, 1, 45, 0, 265] | [70, 500, 1, 25, 0, 345]
2012-07-01 | 2012-09-30 | 92 | 0.26 | 0.33 | 0.35 | 1.74 |  | [170, 250, 1, 45, 0, 265] | [70, 500, 1, 25, 0, 345]
2012-08-01 | 2012-10-31 | 92 | 0.35 | 0.35 | 0.49 | 1.56 |  | [170, 250, 1, 45, 0, 265] | [70, 500, 1, 25, 0, 345]
2012-09-01 | 2012-11-30 | 91 | 0.49 | 0.80 | 0.28 | 1.02 |  | [125, 400, 2, 45, 0, 230] | [50, 150, 4, 50, 0, 290]
2012-10-01 | 2012-12-31 | 92 | 0.80 | 0.93 | 0.25 | 0.75 |  | [170, 250, 1, 45, 0, 265] | [70, 500, 1, 25, 0, 345]
2012-11-01 | 2013-01-31 | 92 | 0.74 | 0.75 | 0.24 | 0.73 |  | [170, 250, 1, 45, 0, 265] | [70, 500, 1, 25, 0, 345]
2012-12-01 | 2013-02-28 | 90 | 0.69 | 0.67 | 0.38 | 0.99 |  | [135, 450, 3, 45, 0, 305] | [110, 200, 1, 10, 0, 255]
2012-01-01 | 2012-06-30 | 181 | 0.54 | 0.48 | 0.62 | 4.80 |  | [145. 400. 5. 45. 0. 265.] | [185. 450. 1. 40. 0. 250.]
2012-02-01 | 2012-07-31 | 181 | 0.44 | 0.45 | 0.60 | 4.10 |  | [140. 400. 2. 45. 0. 310.] | [ 50. 200. 3. 40. 0. 310.]
2012-03-01 | 2012-08-31 | 183 | 0.46 | 0.47 | 0.82 | 5.63 |  | [190. 500. 1. 50. 0. 325.] | [190. 450. 4. 50. 0. 270.]
2012-04-01 | 2012-09-30 | 182 | 0.48 | 0.91 | 0.55 | 7.35 |  | [155. 400. 2. 45. 0. 245.] | [200. 400. 4. 45. 0. 285.]
2012-05-01 | 2012-10-31 | 183 | 0.57 | 0.60 | 1.05 | 8.69 |  | [195. 450. 1. 40. 0. 335.] | [ 50. 250. 1. 45. 0. 320.]
2012-06-01 | 2012-11-30 | 182 | 0.54 | 0.70 | 0.76 | 8.84 |  | [ 90. 400. 1. 45. 0. 310.] | [145. 500. 4. 20. 0. 245.]
2012-07-01 | 2012-12-31 | 183 | 0.80 | 0.77 | 0.70 | 6.67 |  | [ 85. 250. 1. 45. 0. 305.] | [ 70. 500. 1. 25. 0. 345.]
2012-08-01 | 2013-01-31 | 183 | 0.64 | 0.55 | 0.53 | 4.51 |  | [195. 450. 1. 35. 0. 260.] | [160. 400. 1. 35. 0. 245.]
2012-09-01 | 2013-02-28 | 180 | 0.44 | 1.12 | 0.45 | 2.92 |  | [ 65. 450. 4. 45. 0. 215.] | [ 60. 400. 4. 25. 0. 285.]
2012-10-01 | 2013-03-31 | 181 | 0.59 | 1.03 | 1.17 | 4.67 |  | [195. 450. 1. 30. 0. 230.] | [155. 400. 2. 35. 0. 305.]
2012-11-01 | 2013-04-30 | 180 | 0.85 | 0.78 | 1.67 | 8.51 |  | [120. 300. 1. 25. 0. 275.] | [ 60. 250. 4. 50. 0. 255.]
2012-12-01 | 2013-05-31 | 181 | 0.78 | 0.81 | 1.09 | 8.76 |  | [145. 400. 5. 45. 0. 265.] | [ 80. 150. 2. 35. 0. 255.]
2012-02-01 | 2013-01-31 | 366 | 0.71 | 0.74 | 1.13 | 1.32 |  | [155, 450, 4, 45, 0, 350] | [70, 500, 1, 25, 0, 345]
2012-03-01 | 2013-02-28 | 365 | 0.74 | 0.71 | 0.79 | 1.25 |  | [160, 350, 1, 40, 0, 305] | [70, 300, 2, 15, 0, 260]
2012-04-01 | 2013-03-31 | 365 | 0.82 | 0.73 | 0.91 | 1.29 |  | [85, 450, 1, 40, 0, 320] | [120, 400, 4, 25, 0, 270]
2012-05-01 | 2013-04-30 | 365 | 0.71 | 1.10 | 1.27 | 1.35 |  | [185, 450, 1, 45, 0, 255] | [185, 100, 2, 45, 0, 270]
2012-06-01 | 2013-05-31 | 365 | 0.64 | 0.62 | 1.07 | 1.36 |  | [185, 450, 1, 45, 0, 255] | [140, 500, 1, 15, 0, 235]
2012-07-01 | 2013-06-30 | 365 | 0.76 | 0.78 | 1.08 | 1.31 |  | [185, 450, 1, 45, 0, 255] | [85, 400, 5, 25, 0, 240]
2012-08-01 | 2013-07-31 | 365 | 0.86 | 0.92 | 1.11 | 1.25 |  | [110, 400, 1, 45, 0, 260] | [170, 300, 2, 40, 0, 315]
2012-09-01 | 2013-08-31 | 365 | 0.85 | 0.95 | 1.03 | 1.20 |  | [145, 450, 5, 40, 0, 265] | [170, 300, 2, 40, 0, 315]
2012-10-01 | 2013-09-30 | 365 | 0.88 | 0.88 | 1.18 | 1.19 |  | [145, 450, 5, 40, 0, 265] | [140, 300, 1, 45, 0, 260]
2012-11-01 | 2013-10-31 | 365 | 0.82 | 0.85 | 0.87 | 1.19 |  | [145, 450, 5, 40, 0, 265] | [195, 500, 5, 35, 0, 255]
2012-12-01 | 2013-11-30 | 365 | 0.81 | 1.02 | 0.85 | 1.17 |  | [185, 450, 1, 45, 0, 255] | [70, 500, 1, 25, 0, 345]
Bold numbers indicate the best value for a given problem instance. The column
“#obs” indicates the number of missing observations.
“MLP RBF” and “MLP GP” indicate the methods using the MLP with the RBF and GP
model during HPO, respectively.

Figure 4: Daily groundwater level imputation for three months of missing daily
measurement data (November 2012 - January 2013). The missing data have a
linear trend which is captured by imputeTS and the MLP-based method. Figure 5:
Daily Groundwater level imputation for three months of missing daily
measurement data (February 2012 - April 2012). The missing data have a
nonlinear trend and non-interpolating methods like our MLP approach perform
better for these instances. Figure 6: Daily Groundwater level imputation for
three months of missing daily measurement data (March 2012 - May 2012). The
missing data have a nonlinear trend and non-interpolating methods like our MLP
approach perform better for these instances.

In contrast, Figures 5 and 6 show scenarios where the data are behaving
nonlinearly. In Figure 5, the interpolation used by imputeTS failed to capture
the trend and produced an almost flat prediction. Both MLP based methods
captured the nonlinear trend in the missing data better, yet the deviation
over the last $\sim$ 3 weeks was larger. mtsdi made again predictions with
large oscillations around the true data, and for the last 6 weeks of missing
data, it deviated significantly from the truth.

The results for the case where six months of daily observations were missing
are summarized in the middle section of Table 2. Except for two cases (missing
values from July 2012 to December 2012 and from August 2012 to January 2013),
our proposed approach with either RBF or GP made more accurate predictions
than imputeTS and mtsdi in terms of RMSE. Figure 7 is a comparison between all
methods with missing values ranging from November 2012 to April 2013. We can
see that our proposed method captured the high peak better than the other
methods and imputeTS made again a fairly flat prediction. Similarly to the
cases of three months of continuous gaps, mtsdi provided the lowest accuracy.
Note that the increased size of the continuous gap allowed for a larger
probability that seasonality effects will be included in the missing value
ranges.

Figure 7: Daily Groundwater level imputation for six months of missing daily
measurement data (November 2012 - April 2013).

The results for the case where 12 months of daily observations were missing
are summarized in the lower section of Table 2 and in Figure 8. The MLP-based
gap filling method that uses the RBF as surrogate in HPO outperformed both
imputeTS and mtsdi for all instances. Figure 8 shows that the MLP-based gap
filling method captured the nonlinearities in the missing data significantly
better than mtsdi (which introduced strong oscillations and did not capture
the extrema of the missing values) and imputeTS (which failed to capture the
winter data (November 2012-April 2013) and appeared like a piece-wise linear
approximation of the missing data). The MLP-based imputation method showed
stronger deviations from the true data for the last $\sim$1.5 months of gap
filling, which is possibly due to the fact that our method estimates the
missing values sequentially based on the previously estimated values, and thus
it is possible that errors started to accumulate.

Table 2 also show the optimal hyperparameters for the MLP. The results
indicate that the identified optimal MLP architectures depended on the date
range for which the data were missing, and thus the same hyperparameters may
not be optimal for all date ranges and tuning is necessary. The RMSE
performance of a specific architecture was impacted by the data that were
missing. For example, in Table 2, we see that the same optimal hyperparameters
were obtained for the date ranges June 2012-August 2012, July 2012-September
2012, August 2012-October 2012, October 2012-December 2012, and November
2012-January 2013 when using either the RBF or the GP during hyperparameter
tuning, respectively. However, the corresponding RMSEs ranged between 0.26 and
0.8 for RBF, and between 0.33 and 0.93 for the GP. We observe a similar
behavior for the cases of six and twelve months of missing data.

Figure 8: Daily Groundwater level imputation for twelve months of missing
daily measurement data (June 2012 - May 2013). Although all imputation methods
try to capture the nonlinearities, the MLP based imputation agrees best with
the true values.

####  3.3.2 Total water content predictions

Table 3 shows the results of our data imputation for the TWC simulated with
HYDRUS-1D. Shown are the RMSE values in meters between the true and estimated
values obtained with the different gap filling methods for the 12- and
15-month cases. The optimal hyperparameters for the MLP-based approaches are
also shown. Our MLP-based gap filling method achieved lower RMSE values than
ImputeTS for all three tested 12-month date ranges, but it was outperformed by
mtsdi for two test cases. Figure 9 shows the true simulated TWC values and the
values as estimated by the gap filling methods for a 12-month case. We can see
that all methods had difficulties capturing the true data, and our proposed
method and imputeTS did a better job than mtsdi at capturing the first seven
months of missing data, but then both methods generated larger errors during
the remaining five months. In contrast, mtsdi did not capture the trend in the
data of the first seven months well, but the resulting RMSE is smaller overall
regardless.

Table 3: RMSE (in meters) for the different gap filling methods for 12 and 15
months of missing values in daily TWC and optimal hyperparameters ([batch
size, # epochs, # layers, # nodes per layer, dropout rate, # lags]) for both
RBF- and GP-based HPO.

Missing date range (yyyy-mm-dd) | #obs | RMSE |  | Hyperparameters
---|---|---|---|---
From | To |  | MLP RBF | MLP GP | ImputeTS | mtsdi |  | MLP RBF | GP
2013-01-01 | 2013-12-31 | 365 | 0.0084 | 0.0116 | 0.0301 | 0.0155 |  | [50, 500, 1, 50, 0, 45] | [55, 300, 6, 50, 0, 55]
2013-02-01 | 2014-01-31 | 365 | 0.0246 | 0.0233 | 0.0425 | 0.0184 |  | [60, 500, 1, 35, 0, 355] | [70, 500, 1, 25, 0, 345]
2013-03-01 | 2014-02-28 | 365 | 0.0225 | 0.0228 | 0.0364 | 0.0185 |  | [70, 500, 1, 45, 0, 345] | [ 70, 500, 1, 25, 0, 345]
2013-04-01 | 2014-06-30 | 456 | 0.0204 | 0.0210 | 0.0287 | 0.0169 |  | [50, 500, 1, 45, 0, 345] | [70, 500, 1, 25, 0, 345]
2013-05-01 | 2014-07-31 | 457 | 0.0236 | 0.0202 | 0.0280 | 0.0181 |  | [195, 450, 1, 40, 0, 335] | [70, 500, 1, 25, 0, 345]
2013-06-01 | 2014-08-31 | 457 | 0.0233 | 0.0205 | 0.0276 | 0.0187 |  | [195, 450, 1, 40, 0, 335] | [70, 500, 1, 25, 0, 345]
2013-07-01 | 2014-09-30 | 457 | 0.0235 | 0.0209 | 0.0278 | 0.0197 |  | [195, 450, 1, 40, 0, 335] | [70, 500, 1, 25, 0, 345]
2013-08-01 | 2014-10-31 | 457 | 0.0239 | 0.0210 | 0.0273 | 0.0195 |  | [195, 450, 1, 40, 0, 335] | [70, 500, 1, 25, 0, 345]
2013-09-01 | 2014-11-30 | 456 | 0.0053 | 0.0243 | 0.0299 | 0.0178 |  | [50, 500, 1, 50, 0, 30] | [ 70, 500, 1, 25, 0, 345]
2013-10-01 | 2014-12-31 | 457 | 0.0276 | 0.0282 | 0.0383 | 0.0235 |  | [195, 450, 1, 40, 0, 335] | [70, 500, 1, 25, 0, 345]
2013-11-01 | 2015-01-31 | 457 | 0.0291 | 0.0266 | 0.0472 | 0.0229 |  | [195, 450, 1, 40, 0, 335] | [70, 500, 1, 25, 0, 345]
2013-12-01 | 2015-02-28 | 455 | 0.0141 | 0.0203 | 0.0308 | 0.0228 |  | [50, 450, 1, 45, 0, 120] | [70, 500, 1, 25, 0, 345]
Bold numbers indicate the best value for a given problem instance. The column
“#obs” indicates the number of missing observations.
“MLP RBF” and “MLP GP” indicate the methods using the MLP with the RBF and GP
model during HPO, respectively.

Table 3 also shows the results for the case of 15 months of missing data. The
MLP-based gap filling method outperformed imputeTS for all instances, but it
was again outperformed by mtsdi for most cases (mtsdi had a smaller RMSE).
Figure 10 shows the results for the gap September 2013-November 2014 (a case
in which the MLP based method performed better than the other methods). We can
see that the MLP that used the RBF in the HPO led to predictions that almost
perfectly matched the true data. mtsdi appeared to be able to capture some of
the trends in the data, in particular the large peaks with higher frequency
oscillations, but it failed to approximate the “less rugged” and lower values
well. We can also observe from Figure 10 how important the choice of
hyperparameters is for the MLP. The MLP using the GP during HPO did not
perform nearly as well as the MLP that used the RBF during HPO. The main
differences between both architectures were in the batch size (50 vs. 70), the
number of nodes per layer (50 vs. 25), and the number of lags (30 vs. 345). On
the other hand, in Figure 9, the GP-based and RBF-based MLP predictions were
approximately equal, and for this case the optimal architectures differed only
with respect to the number of nodes per layer (45 vs. 25), with the larger
architecture (RBF solution) being marginally better than the GP based
solution. Thus, we can recognize the importance of HPO when using DL models
and the sensitivity of the MLP’s performance with regard to the architecture.

Figure 9: Comparison of gap filling methods for daily TWC with missing values
between March 2013 and February 2014. Figure 10: Comparison of gap filling
methods for daily TWC with missing values between September 2013 and November
2014.

Figure 11 shows the imputation results with 15 months of missing values from
April 2013 to June 2014. Similar to Figure 9, mtsdi provided the smallest RMSE
but our proposed MLP based methods provided better estimates at the beginning
while mtsdi overestimated the values by making an almost constant value
prediction.

Figure 11: Comparison of gap filling methods for daily TWC with missing values
between April 2013 and June 2014.

####  3.3.3 FLUXNET predictions

Table 4 shows the results for the FLUXNET2015 dataset. Because the number of
available observations was large, we stress-tested the performance of the gap
filling methods for very large continuous data gaps, using between 2,159 (3
months) and 35,063 (36 months) missing observations. The results show that the
MLP-based gap filling method attained better performance (in terms of RMSE)
than imputeTS and mtsdi. The top panel of Figure 12 shows the true and the
estimated values of NEE_VUT_USTAR50 for the case of 8,759 missing values. The
lower panels show scatter plots of the true values versus the predicted
values. If the methods were able to exactly repredict the missing values, all
points would lie on the diagonal. As can be seen in Figure 12, imputeTS failed
to capture the variability in the data and it simply connected the last
observation (last data point before the gap) with the first observation after
the gap. The MLP-based methods on the other hand made predictions that were
close to the true values (points lie close to the diagonal). The point cloud
obtained with mtsdi shows a larger variability around the diagonal, and thus
for this use case mtsdi made less reliable predictions. Figure 13 plots the
results for the case when 17,543 observations are missing. Again, imputeTS
failed to capture the trend in the data completely. mtsdi captured the
seasonal trend but it failed to accurately predict the large and small values.
Our proposed MLP-based imputation methods estimated the high and low values
significantly better. This is also reflected in the scatter plots which show
that the point clouds of our proposed method lie closer to the diagonal than
those of imputeTS and mtsdi.

Comparing the optimal hyperparameters found by the RBF and the GP methods,
respectively, we can see that very different MLP architectures can lead to
very similar RMSEs. For example, for the case where 2,159 observations were
missing, the difference in RMSEs was only 0.02, but the GP used one extra
hidden layer and a significantly larger number of lags. On the other hand, the
RBF-based HPO method used more epochs. Similarly, for the case where 17,543
observations were missing, the RBF-based optimizer led to a significantly
larger architecture, but it improved the RMSE by only 0.02. Thus, there is a
trade-off between architecture size (and thus training time) and RMSE, and
perhaps for large datasets, accepting a slightly larger RMSE in exchange for a
simpler and faster to train network is preferable.

Table 4: RMSE for different gap filling methods for hourly FLUXNET dataset and
hyperparameters ([batch size, # epochs, # layers, # nodes per layer, dropout
rate, # lags]) for both RBF- and GP-based HPO.

Missing Range | # obs | RMSE |  | Hyperparameters
---|---|---|---|---
From | To | RBF | GP | ImputeTS | mtsdi |  | RBF | GP
2007-01-01 | 2007-03-31 | 2,159 | 0.98 | 0.96 | 1.86 | 4.83 |  | [190, 450, 2, 45, 0, 50] | [180, 100, 3, 40, 0, 240]
2007-01-01 | 2007-12-31 | 8,759 | 2.72 | 2.68 | 8.14 | 5.17 |  | [ 55, 400, 4, 40, 0, 30] | [150, 200, 1, 35, 0, 355]
2007-01-01 | 2008-12-31 | 17,543 | 2.86 | 2.88 | 8.20 | 5.14 |  | [ 60, 350, 5, 45, 0, 90] | [55, 200, 2, 35, 0, 185]
2007-01-01 | 2010-12-31 | 35,063 | 2.75 | 2.94 | 8.15 | 5.24 |  | [ 60, 350, 5, 45, 0, 90] | [70, 500, 1, 25, 0, 345]
Bold numbers indicate the best value for a given problem instance. The column
“#obs” indicates the number of missing observations.
“MLP RBF” and “MLP GP” indicate the methods using the MLP with the RBF and GP
model during HPO, respectively.

Figure 12: Hourly NEE_VUT_USTAR50 (FLUXNET) imputation comparison: missing
value range is January 1, 2007 to December 31, 2007. The top panel shows the
true and the predicted values. The bottom panels show scatter plots for each
method (ideally, prediction and true data are identical and lie on the
diagonal)

Figure 13: Hourly NEE_VUT_USTAR50 imputation comparison: missing value range
2007-01-01–2008-12-31

####  3.3.4 Discussion

The three gap filling case studies showed that there was no one single method
that performed best for all problems and all instances. We found that for the
two examples with real observation data (as opposed to simulation-created
data), our proposed gap filling method tended to outperform both ImputeTS and
mtsdi in terms of RMSE, in particular when the gaps were large. As described
in the previous subsections, ImputeTS performed well when the missing data can
be approximated by a somewhat (piecewise) linear function. However, for the
FLUXNET2015 data, ImputeTS failed, which may be due to the amount of missing
data, or perhaps due to the high frequency oscillations in the data. mtsdi
performed better than our method for the simulated total water content data in
terms of RMSE, although it did not seem to capture much of the details of the
lowest values or the dynamics in the data. Figures 9 and 11 show that mtsdi
seemed to capture the large faster oscillating peaks quite well, but for the
slower changes of TWC, it appeared to use a mean value instead of capturing
the slower variabilities in the data.

Although our proposed imputation approach performed generally well and often
outperformed the other methods, its computing time requirement was
significantly larger than that of imputeTS and mtsdi. This is due to the
hyperparameter optimization trying many different sets of hyperparameters and
each training is computationally nontrivial. On the other hand, if the values
to be gap filled behave nonlinearly, and if high prediction accuracy is
desired, we recommend using the DL based approach due to the limitations of
the other methods in these cases. In practice, one would not conduct a
detailed study as presented here with artificially created gaps, but rather
optimize and train the MLP once for the time series to be gap filled. One way
to speed up the DL approach is to either parallelize the repetitions of the
training if compute resources are available or to reduce the number of
repetitions. However, the latter option may lead to models that are less
reliable in terms of predictive performance. In the FLUXNET test cases, we
used only one training repetition (instead of five as done with the two other
use cases). The MLP based gap filling method always performed significantly
better than the other two methods, indicating that it was not simply a “lucky”
draw.

Many environmental datasets contain functional relationships between data
features. However, many time series imputation methods (including the methods
ImputeTS and mtsdi discussed in this paper) do not take advantage of this
additional information and they use only statistical information contained in
the time series data. In other words, they are purely data-driven approaches.
Unlike these approaches, our proposed DL method can be extended to include any
known physical relationships and constraints in addition to the statistical
information. For example, in the HYDRUS-1D case, we know that there is a mass
conservation constraint present between the features and it is possible to
include this constraint, e.g., by adding a penalty term to the objective
function or by including hard constraints during model training.

In our study, we compared the different methods based on their RMSE
performance. Although the RMSE is a widely used measure, it does not reflect
the individual errors well. For example, in the daily TWC imputation
comparison shown in Figure 9, one could argue that the proposed MLP based
methods produce fairly accurate and better predictions than mtsdi for 50% of
the missing values, but they accumulate larger errors for the remaining
values. On the other hand, mtsdi does not match most of the data well, but its
almost flat prediction accrues overall a lower RMSE. Thus, it is possible that
optimizing for another performance measure may lead to better outcomes for the
MLP methods.

##  4 Summary and Conclusions

In this paper, we proposed a time series missing value imputation method that
uses DL models. Our focus was on time series imputation problems where a long
continuous gap was present in one variable, and where all other supporting
variables are fully observed. Our proposed method uses an MLP whose
architecture is tuned with a derivative-free surrogate model based
optimization method. To this end, we modified the hyperparameter optimization
approach proposed in muller2019surrogate  in order to facilitate the
imputation of missing values in time series.

After training the MLP, we imputed the missing values sequentially from the
first missing value to the last missing value, basing predictions on previous
predictions as we go along. Our sequential estimation of missing values
allowed us to fill long-term gaps. We chose MLPs as DL models because they
have previously shown good performance for time series data and because the
time needed to optimize the model is significantly shorter than for other
types of DL models (see muller2019surrogate ). Generally, however, any
feedforward artificial neural network can be used within our proposed method
such as convolutional neural networks fukushima1980neocognitron ;
lecun1989backpropagation ; lecun1998gradient .

We performed numerical experiments with three different test cases in which we
created artificial gaps which allowed us to assess the quality of the gap
filled values. We used two test cases with observed values and one test case
with simulated data. The numerical results showed that for many problem
instances, our proposed approach was able to fill the gaps with highly
accurate values and it was able to capture trends and dynamics in the data
better than other data-driven methods such as ImputeTS and mtsdi. Using the
RBF as surrogate model during HPO appeared to yield more stable performance
than using the GP. Therefore, we recommend using our proposed approach with
the RBF surrogate during HPO.

The numerical results also showed that the optimal MLP hyperparameters
depended on which data were missing (gap size and location in the time
series). Therefore, hyperparameter optimization is essential before using any
ML model to impute the missing values.

The results of this paper can lay the foundation for studying new approaches
for filling long gaps in time series. Although the proposed method was able to
fill large gaps with seasonal trends well, including capturing large and small
values, it sometimes failed to correctly estimate the missing values towards
the end of the data gap, indicating a potential for accumulation of errors as
we inform predictions with previously predicted values. A “backward imputation
approach” might be helpful to alleviate this drawback. In this approach, we
could modify the MLP inputs such that we start imputing data from the end of
the gap backwards in time. Combining both forward and backward imputation may
allow us to reap the benefits of both.

Our study focused on time series that have one large gap. However, we expect
the method to perform reasonably well when we have multiple large gaps in the
time series. Since the inputs to the MLP are designed such that they are
independent and the HPO is general, we expect that filling multiple gaps with
this method is straight-forward. For time series with a mix of large and small
gaps, a combination of our DL model based method and other methods that are
aimed at filling small gaps can be used. Incorporation of physics constraints
in the DL model predictions has the potential to improve prediction accuracy
of environmental data analysis and predictions. Moreover, if values are
missing in several variables of the multivariate time series, it may be
possible to modify the proposed approach for this scenario. This is a future
research topic.

This paper studied the time-series imputation with the MLP which is the most
basic feedforward neural network. It would also be interesting to study the
proposed time-series imputation approach with other advanced feedforward
neural networks such as convolutional neural network lecun1998gradient  and
AR-NET triebe2019ar . Recently, a state-of-art time-series forecasting
software NeuralProphet triebe2021neuralprophet  combines the AR-NET and the
Prophet software taylor2018forecasting .

## Data availability

Three datasets used in this study have been deposited in the GitHub public
repository: https://github.com/JanghoPark-
LBL/Dataset_Missing_Value_Imputation_with_Deep_Neural_Networks.

## Acknowledgements

This work was supported by Laboratory Directed Research and Development (LDRD)
funding from Berkeley Lab, provided by the Director, Office of Science, of the
U.S. Department of Energy under Contract No. DE-AC02-05CH11231. This research
used resources of the National Energy Research Scientific Computing Center
(NERSC), a U.S. Department of Energy Office of Science User Facility operated
under Contract No. DE-AC02-05CH11231.

## References

  * (1) García-Laencina, P.J., Sancho-Gómez, J.-L., Figueiras-Vidal, A.R.: Pattern classification with missing data: a review. Neural Computing and Applications 19(2), 263–282 (2010)
  * (2) Yozgatligil, C., Aslan, S., Iyigun, C., Batmaz, I.: Comparison of missing value imputation methods in time series: the case of Turkish meteorological data. Theoretical and Applied Climatology 112(1-2), 143–167 (2013)
  * (3) Kalteh, A.M., Hjorth, P.: Imputation of missing values in a precipitation–runoff process database. Hydrology Research 40(4), 420–432 (2009)
  * (4) Aissia, M.-A.B., Chebana, F., Ouarda, T.B.: Multivariate missing data in hydrology–review and applications. Advances in Water Resources 110, 299–309 (2017)
  * (5) Gu, C., Anderson, W., Maggi, F.: Riparian biogeochemical hot moments induced by stream fluctuations. Water Resources Research 48(9) (2012)
  * (6) Arora, B., Wainwright, H.M., Dwivedi, D., Vaughn, L.J., Curtis, J.B., Torn, M.S., Dafflon, B., Hubbard, S.S.: Evaluating temporal controls on greenhouse gas (ghg) fluxes in an arctic tundra environment: An entropy-based approach. Science of the total environment 649, 284–299 (2019)
  * (7) Thi-Thu-Hong Phan, B., Caillault, E.P., Bigand, A.: edtwbi: Effective imputation method for univariate time series. In: Advanced Computational Methods for Knowledge Engineering: Proceedings of the 6th International Conference on Computer Science, Applied Mathematics and Applications, ICCSAMA 2019, vol. 1121, p. 121 (2019). Springer Nature
  * (8) Moritz, S., Sardá, A., Bartz-Beielstein, T., Zaefferer, M., Stork, J.: Comparison of different methods for univariate time series imputation in R. arXiv preprint arXiv:1510.03924 (2015)
  * (9) Moritz, S., Bartz-Beielstein, T.: imputeTS: time series missing value imputation in R. The R Journal 9(1), 207–218 (2017)
  * (10) Phan, T.-T.-H.: Machine learning for univariate time series imputation. In: 2020 International Conference on Multimedia Analysis and Pattern Recognition (MAPR), pp. 1–6 (2020)
  * (11) Batista, G.E.A.P.A., Monard, M.C.: A study of K-nearest neighbour as an imputation method. Hybrid Intelligent Systems 87(251-260), 48 (2002)
  * (12) Che, Z., Purushotham, S., Cho, K., Sontag, D., Liu, Y.: Recurrent neural networks for multivariate time series with missing values. Scientific Reports 8(1), 1–12 (2018)
  * (13) Cao, W., Wang, D., Li, J., Zhou, H., Li, L., Li, Y.: Brits: Bidirectional recurrent imputation for time series. In: NeurIPS (2018)
  * (14) Luo, Y., Cai, X., Zhang, Y., Xu, J., et al.: Multivariate time series imputation with generative adversarial networks. In: Advances in Neural Information Processing Systems, pp. 1596–1607 (2018)
  * (15) Zhang, Y., Zhou, B., Cai, X., Guo, W., Ding, X., Yuan, X.: Missing value imputation in multivariate time series with end-to-end generative adversarial networks. Information Sciences 551, 67–82 (2021)
  * (16) Lingras, P., Zhong, M., Sharma, S.: Evolutionary regression and neural imputations of missing values. In: Soft Computing Applications in Industry, pp. 151–163 (2008)
  * (17) Zhong, M., Sharma, S., Lingras, P.: Rationalizing reliable imputation durations of genetically designed time delay neural network and locally weighted regression models. Transportation Planning and Technology 30(6), 609–626 (2007)
  * (18) Mahabbati, A., Beringer, J., Leopold, M., McHugh, I., Cleverly, J., Isaac, P., Izady, A.: A comparison of gap-filling algorithms for eddy covariance fluxes and their drivers. Geoscientific Instrumentation, Methods and Data Systems Discussions, 1–31 (2020)
  * (19) Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cournapeau, D., Brucher, M., Perrot, M., Duchesnay, E.: Scikit-learn: Machine learning in Python. Journal of Machine Learning Research 12, 2825–2830 (2011)
  * (20) Kim, Y., Johnson, M.S., Knox, S.H., Black, T.A., Dalmagro, H.J., Kang, M., Kim, J., Baldocchi, D.: Gap-filling approaches for eddy covariance methane fluxes: A comparison of three machine learning algorithms and a traditional method with principal component analysis. Global Change Biology 26(3), 1499–1518 (2020)
  * (21) Waibel, A., Hanazawa, T., Hinton, G., Shikano, K., Lang, K.J.: Phoneme recognition using time-delay neural networks. IEEE Transactions on Acoustics, Speech, and Signal Processing 37(3), 328–339 (1989)
  * (22) Müller, J., Park, J., Sahu, R., Varadharajan, C., Arora, B., Faybishenko, B., Agarwal, D.: Surrogate optimization of deep neural networks for groundwater predictions. Journal of Global Optimization, 10–100710898020009120 (2020)
  * (23) Junger, W., Ponce de Leon, A.: Package ’mtsdi’
  * (24) Nair, V., Hinton, G.E.: Rectified linear units improve restricted boltzmann machines. ICML’10, pp. 807–814. Omnipress, Madison, WI, USA (2010)
  * (25) Kingma, D.P., Ba, J.: Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980 (2014)
  * (26) Rumelhart, D.E., Hinton, G.E., Williams, R.J.: Learning representations by back-propagating errors. Nature 323(6088), 533–536 (1986)
  * (27) Srivastava, N., Hinton, G., Krizhevsky, A., Sutskever, I., Salakhutdinov, R.: Dropout: a simple way to prevent neural networks from overfitting. The Journal of Machine Learning Research 15(1), 1929–1958 (2014)
  * (28) Goodfellow, I., Bengio, Y., Courville, A., Bengio, Y.: Deep Learning vol. 1. MIT press Cambridge, ??? (2016)
  * (29) Feurer, M., Hutter, F.: Hyperparameter optimization. In: Automated Machine Learning, pp. 3–33. Springer, ??? (2019)
  * (30) Bergstra, J., Bengio, Y.: Random search for hyper-parameter optimization. Journal of Machine Learning Research 13(2) (2012)
  * (31) Larochelle, H., Erhan, D., Courville, A., Bergstra, J., Bengio, Y.: An empirical evaluation of deep architectures on problems with many factors of variation. In: Proceedings of the 24th International Conference on Machine Learning, pp. 473–480 (2007)
  * (32) Snoek, J., Larochelle, H., Adams, R.P.: Practical Bayesian optimization of machine learning algorithms. arXiv preprint arXiv:1206.2944 (2012)
  * (33) Xie, L., Yuille, A.: Genetic CNN. In: Proceedings of the IEEE International Conference on Computer Vision, pp. 1379–1388 (2017)
  * (34) Sahu, R.K., Müller, J., Park, J., Varadharajan, C., Arora, B., Faybishenko, B., Agarwal, D.: Impact of input feature selection on groundwater level prediction from a multi-layer perceptron neural network. Frontiers in Water 2, 46 (2020)
  * (35) Šimŭnek, J., Šejna, M., Saito, H., Sakai, M., van Genuchten, M.T.: The HYDRUS-1D Software Package for Simulating the One-Dimensional Movement of Water, Heat, and Multiple Solutes in Variably-Saturated Media, Version 4.17, edn. Riverside, California, (2013). Riverside, California
  * (36) Šimŭnek, J., van Genuchten, M.T.: Modeling nonequilibrium flow and transport processes using hydrus. Vadose Zone Journal 7(2), 782–797 (2008)
  * (37) Arora, B., Mohanty, B.P., McGuire, J.T.: An integrated Markov chain Monte Carlo algorithm for upscaling hydrological and geochemical parameters from column to field scale. Science of the Total Environment 512, 428–443 (2015)
  * (38) Baek, S., Ligaray, M., Pachepsky, Y., Chun, J.A., Yoon, K.-S., Park, Y., Cho, K.H.: Assessment of a green roof practice using the coupled SWMM and HYDRUS models. Journal of Environmental Management 261, 109920 (2020)
  * (39) Pastorello, G., Trotta, C., Canfora, E., Chu, H., Christianson, D., Cheah, Y.-W., Poindexter, C., Chen, J., Elbashandy, A., Humphrey, M., et al.: The FLUXNET2015 dataset and the ONEFlux processing pipeline for eddy covariance data. Scientific Data 7(1), 1–27 (2020)
  * (40) United States Department of Energy, T.O.o.S.: FLUXNET2015 Dataset. https://fluxnet.org/data/fluxnet2015-dataset/. Last accessed: 2021-1-1
  * (41) Novick, K., Phillips, R.: FLUXNET2015 US-MMS Morgan Monroe State Forest. FLUXNET; Indiana University, https://doi.org/10.18140/flx/1440083 (2016)
  * (42) Paszke, A., Gross, S., Chintala, S., Chanan, G., Yang, E., DeVito, Z., Lin, Z., Desmaison, A., Antiga, L., Lerer, A.: Automatic differentiation in PyTorch. In: 31st Conference on Neural Information Processing Systems (NIPS 2017), Long Beach, CA, USA (2017)
  * (43) Fukushima, K.: Neocognitron: A self-organizing neural network model for a mechanism of pattern recognition unaffected by shift in position. Biological Cybernetics 36(4), 193–202 (1980)
  * (44) LeCun, Y., Boser, B., Denker, J.S., Henderson, D., Howard, R.E., Hubbard, W., Jackel, L.D.: Backpropagation applied to handwritten zip code recognition. Neural Computation 1(4), 541–551 (1989)
  * (45) LeCun, Y., Bottou, L., Bengio, Y., Haffner, P.: Gradient-based learning applied to document recognition. Proceedings of the IEEE 86(11), 2278–2324 (1998)
  * (46) Triebe, O., Laptev, N., Rajagopal, R.: Ar-net: A simple auto-regressive neural network for time-series. arXiv preprint arXiv:1911.12436 (2019)
  * (47) Triebe, O., Hewamalage, H., Pilyugina, P., Laptev, N., Bergmeir, C., Rajagopal, R.: Neuralprophet: Explainable forecasting at scale. arXiv preprint arXiv:2111.15397 (2021)
  * (48) Taylor, S.J., Letham, B.: Forecasting at scale. The American Statistician 72(1), 37–45 (2018)


# A Tutorial on the Pretrain-Finetune Paradigm for Natural Language
Processing111Our replication package can be accessed at
https://doi.org/10.7910/DVN/QTT84C.

Yu Wang
Institute for Advanced Study in Social Sciences
Fudan University
<EMAIL_ADDRESS>
###### Abstract

The pretrain-finetune paradigm represents a transformative approach in natural
language processing (NLP). This paradigm distinguishes itself through the use
of large pretrained language models, demonstrating remarkable efficiency in
finetuning tasks, even with limited training data. This efficiency is
especially beneficial for research in social sciences, where the number of
annotated samples is often quite limited. Our tutorial offers a comprehensive
introduction to the pretrain-finetune paradigm. We first delve into the
fundamental concepts of pretraining and finetuning, followed by practical
exercises using real-world applications. We demonstrate the application of the
paradigm across various tasks, including multi-class classification and
regression. Emphasizing its efficacy and user-friendliness, the tutorial aims
to encourage broader adoption of this paradigm. To this end, we have provided
open access to all our code and datasets. The tutorial is particularly
valuable for quantitative researchers in psychology, offering them an
insightful guide into this innovative approach.

Keywords: Large language models, pretrain-finetune, natural language
processing

##  1 Introduction

The rise of pretrain-finetune paradigm has greatly reshaped the landscape of
natural language processing (NLP) (E. Hu ., 2022). This new paradigm,
characterized by the application of large pretrained language models, most
notably BERT (Devlin ., 2019) and RoBERTa (Liu ., 2019), and high efficacy on
finetuned tasks even in the face of relatively few training samples, is now
being widely applied in social sciences (Zhang ., 2021; Wang, 20231, 20232).
While earlier works have recommended a minimum of 3,000 samples for NLP tasks
using the bag-of-words approach (Wang ., 2022), recent research applying the
pretrain-finetune paradigm has shown that finetuned large models with just a
few hundred labeled samples could yield competitive performance (Wang, 20232).

Figure 1: High-level illustration of the pretraining (left) and finetuning
(right) workflows.

In the current tutorial, we aim to provide (1) an overview of the pretrain-
finetune paradigm and (2) illustrative applications of the pretrain-finetune
paradigm to research questions in social sciences (Figure 1). We first provide
an introduction to the key concepts in pretraining and finetuning, such as
tokenization and encoding. We then use some of the most recent applications of
the pretrain-finetune paradigm in social sciences to illustrate how to use
this paradigm to advance quantitative research in our discipline.

##  2 Pretrain

Pretraining is the process of training a model with unlabeled raw data with no
particular downstream tasks. Some of the most widely used pretrained models
include BERT and RoBERTa. These large language models usually contain hundreds
of millions of parameters. Pretraining these models from scratch requires
access to large amounts of raw data and specialized hardware like graphics
processing units (GPUs). The pretraining process consists of the following
steps: (1) tokenization, (2) encoding, and (3) pretraining tasks, such as
masked language modeling, next sentence prediction (Devlin ., 2019) and
sentence-order prediction (Lan ., 2020).

###  2.1 Tokenization

Unlike earlier methods such as bag of words (Wang ., 2022) or word embeddings
(Mikolov ., 2013), large langauge models mostly use subwords as a token, such
as wordpieces in BERT (Devlin ., 2019) and Byte-Pair Encoding (BPE) in RoBERTa
(Liu ., 2019) and GPT-3 (Brown ., 2020). To illustrate how words are broken
into subwords, let’s use a text snippet from the Journal’s description and
tokenize it using the BPE tokenizer used in GPT-3.5 and GPT-4.222The snippet
is taken from https://journals.sagepub.com/description/AMP. The tokenizer is
available at https://platform.openai.com/tokenizer and was accessed on January
6th 2024..

    AMPPS encourages integration of methodological and analytical questions and
    brings the latest methodological advances to non-methodology experts across
    all areas of the field.

The word “AMPPS” is split into “AM”, “PP” and “S”. The word “methodological”
is split into “method”, “ological”. And lastly the word “non-methodology” is
split into “non”, “-method”, and “ology”. All other words, including
punctuation marks, remain individual units without further splitting. These
units are then turned into tokens represented with non-negative integers.

    AM PP S encourages integration of method ological and analytical questions
    and brings the latest method ological advances to non -method ology
    experts across all areas of the field .

These units are then turned into tokens represented with integers.

    [1428, 4505, 50, 37167, 18052, 315, 1749, 5848, 323, 44064, 4860, 323, 12716,
    279, 5652, 1749, 5848, 31003, 311, 2536, 51397, 2508, 11909, 4028, 682, 5789,
    315, 279, 2115, 627]

###  2.2 Encoding

Once we have tokenized the input text into tokens, we then retrieve the
corresponding embeddings for these tokens, where the token id, a non-negative
integer, would serve the retrieval key. These token embeddings are vector
representations. For the BERT-large model, each such vector contains 1024
float numbers. These sets of vectors are then fed into the encoder
layers.333BERT models use encoders and are the focus of this tutorial. Note
that GPT models use decoders.

Figure 2: Illustration of the BERT model’s self-attention mechanism.

The key component of the encoders is self-attention (Vaswani ., 2017).
Intuitively, it computes how much attention each token pays to the other
tokens and generates the new representation for a token by calculating the
weighted average of all the tokens where the weight is based on self-
attention.444For a similar illustration, readers can refer to Kjell . (2023).
For a more in-depth illustration, please refer to the original paper (Vaswani
., 2017) and the tutorial on transformers at
http://jalammar.github.io/illustrated-transformer/. Such encoders are then put
on top of each other in a sequential manner. For example, in BERT-large models
there are 24 encoder layers. In BERT-base models there are 12 encoder layers.

###  2.3 Pretraining tasks

Now that we have a good understanding of tokenization and encoding, let’s
proceed to discuss how these parameters in these layers are trained using
pretraining tasks. The goal is to train these large language models from
scratch and embue them with world knowledge that exists in the training corpus
(Longpre ., 2020). Commonly used pretraining tasks include masked language
modeling, next sentence prediction (Devlin ., 2019) and sentence order
prediction (Lan ., 2020).555Later research has shown that the next sentence
prediction does not contribute much to the pretraining process (Liu ., 2019).

Figure 3: Illustration of the BERT model’s masked language modeling.

Let’s see an example of how masked language modeling works. Suppose we have
the following input text: “I enjoy trying new activities”. During training,
some of the tokens will be randomly masked and the language model is tasked
with predicting these masked tokens. In Figure 3, the token “activities” is
masked. The language model then is tasked with predicting this masked token.
All tokens in the vocabulary are eligible candidates. Suppose the model
predicts “activities” with a probability of 0.2. If we use cross entropy loss,
then the loss term for this prediction is -log(0.2), which is its negative
log-likelihood.666For details, please see https://github.com/google-
research/bert/blob/master/run_pretraining.py#L303. Note that the higher
probability the model assigns to the corrected token, the lower the loss. The
language model is trained on this task over the entire corpus for a few times
until the loss in prediction stops decreasing.

##  3 Finetuning

Finetuning is the process of training an already-pretrained model on a
downstream task. It is characterized by (1) adding a small set new parameters
to accommodate the new tasks and (2) a relatively small learning rate given
the fact that the vast majority of the parameters are already reasonably
trained during the pretraining stage.777In terms of training procedures, it is
the same as training a supervised model from scratch. Readers could refer to a
recent tutorial on how to train supervised models by Pargent . (2023).
Downstream tasks at this stage can be broadly grouped into classification and
regression. Examples of classification tasks include sentiment analysis and
topic classification; examples of regression include fatality prediction. This
is the step where quantitative psychologists need to apply a pretrained model
to a particular task, such as depression detection and personality
classification. For each specific task, researchers need to specify the
model’s input, targeted output (including the number of labels), and whether
this is a classification task or a regression task.

Figure 4: Illustration of the finetuning process. Same as with pretraining,
the model takes in a text input and generates an encoding for each token. What
differentiates finetuning from pretraining is that it takes a particular token
representation to make task-specific predictions.

In Figure 4, we illustrate how a k-class classification works. From the
encoder layers, we retrieve a 1 by N vector for each input sample, where N is
768 in the case of BERT-base and 1024 in the case of BERT-large. Then we
project this vector to 1 by K using an N by K matrix. This layer is often
referred to as the fully connected layer. What follows next depends on whether
the task is classification or regression. In the case of classification, where
K is equal to or greater than 2, we apply a softmax function such that each of
the K elements represents the probability of that corresponding class being
the predicted class:888For implementation details, refer to
https://pytorch.org/docs/stable/generated/torch.nn.Softmax.html.

$softmax(x_{i})=\frac{exp(x_{i})}{\sum_{j}(exp(x_{j}))}$

In the case of regression, where K is 1, we can apply the loss function
directly to the 1 by K vector (now 1 by 1). Mean squared error is a commonly
used loss function for regression tasks.999For implementation details, refer
to https://pytorch.org/docs/stable/generated/torch.nn.MSELoss.html.

##  4 Practical Exercises

In this section, we provide two practical exercises to illustrate how
researchers can finetune a large language model with relative ease to achieve
state-of-the-art results in classification and regression,
respectively.101010Some of the original results on classification and
regression have been published at Political Analysis. All our exercises are
written in Python in the format of Jupyter notebooks so that readers can
follow along in an interactive manner.111111https://jupyter.org. For easy
reproducibility and access to computation resources, all our computation is
done on Google Colab.121212https://colab.research.google.com.

###  4.1 Multi-class classification

Figure 5:  Given an input text, we categorize it into one of eight mutually
exclusive classes. In the provided example, the text is classified under the
second topic, External Relations.

When researchers need to classify text inputs into more than two classes, this
is what we call multi-class classification. In the following case study, we
classify a given text into one of eight topics. The texts are the speech
transcripts from the New Zealand Parliament from 1987 to 2002 (Osnabrügge .,
2021). In Figure 5, we illustrate the data input and list all the eight
classes. We use using 2,915 annotated speech transcripts for training, 625 for
validation, and another 625 for testing.

We finetune a RoBERTa-base model (Liu ., 2019) for topic classification.
RoBERTa-base has 12 layers of transformers and 125 million parameters in
total. On top of its 12 layers of transformers, we add a classification layer
for 8-topic classification. We use cross entropy as the loss function. We
finetune the RoBERTa-base model for 20 epochs with a learning rate of 2e-5, a
batch size of 16, and an input sequence length of 512 on an A100 GPU. We use
the validation set’s cross entropy loss to select the best epoch and the
optimal checkpoint. We then use the optimal checkpoint to make inferences on
the test set with a batch size of 64 (see Listing 1). As our baseline, we use
results from Osnabrügge . (2021), where the authors use 115,410 related but
out-of-domain policy statements from Australia, Canada, Ireland, New Zealand,
the United Kingdom, and the United States to train a regularized multinomial
logistic regression model. For easy comparison, we use the same evaluation
metrics as in Osnabrügge . (2021).

⬇

1

2 training_args = TrainingArguments(

3 output_dir="./results", # output directory

4 num_train_epochs=epochs, # total number of training epochs

5 per_device_train_batch_size=16, # batch size for training

6 per_device_eval_batch_size=64, # batch size for evaluation

7 learning_rate = 2e-5, # learning rate

8 save_strategy= "epoch", # save checkpoints after each epoch

9 evaluation_strategy="epoch",

10 load_best_model_at_end= True,

11 )

12

13 def model_init():

14 return RobertaForSequenceClassification.from_pretrained("roberta-base",
num_labels=8) # number of labels

Listing 1: Specify the number of labels to be 8 for the multi-class
classification and set the hyperparameters such as batch size and learning
rate.

Across all metrics, finetuning the RoBERTa model with 2,915 New Zealand
parliamentary speeches substantially outperforms the cross-domain topic
classifier by Osnabrügge . (2021), which is trained using 115,420 annotated
policy statements (Table 1). For example, for top-1 accuracy the finetuned
RoBERTa model outperforms the cross-domain baseline by 26%. For top-3 accuracy
the finetuned RoBERTa model outperforms the cross-domain baseline by 10%.

Table 1: Finetuning a RoBERTa-base model with 2,915 labeled in-domain samples can outperform the cross-domain regularized multinomial logistic regression model, trained with 115,420 out-of-domain samples, in the 8-topic classification task by a large margin. Cross-domain classifiers are from Osnabrügge . (2021). Test set is the same for both models. Mean of three random runs is reported, with standard deviation in the brackets. Better results are in bold. Metrics | 8 topics
---|---
Cross-domain | Finetuning LM
Top-1 accuracy | 0.512 (0.007) | 0.643 (0.007)
Top-3 accuracy | 0.821 (0.001) | 0.899 (0.001)
Top-5 accuracy | 0.917 (0.007) | 0.968 (0.007)
Balanced accuracy | 0.465 (0.004) | 0.592 (0.004)
F1 macro | 0.456 (0.011) | 0.584 (0.011)

###  4.2 Regression

In regression tasks, researchers use an input text to predict a numerical
value. In the following exercise, we predict the natural logarithm of
fatalities on a country month level using the CrisisWatch text (Häffner .,
2023). For a visual representation of the input and output, please refer to
Figure 6.

Figure 6: Given an input text, we predict a numerical value. In the example
above, given an expert-written International Crisis Group (ICG) CrisisWatch
report on Afghanistan for October, 2006, the model is expected to predict 6.34
as the natural logarithm of fatalities.

In Häffner . (2023), for the task of predicting the natural logarithm of
fatalities on a country month level using the CrisisWatch text, the authors
use random forest and XGBoost.131313For an introduction to these tree-based
models, readers can refer to Pargent . (2023) and Wang (2019). The training
set covers the period from 2003 to first half of 2020 and has 21,924 samples,
the validation set covers the second half of 2020 and has 648 samples, and the
test set covers the year 2021 and has 1296 samples. These two models achieve
an R2 of 0.64 and 0.63 respectively and an MSE of 1.59 and 1.60 respectively.
We use the same training set, validation set and testing set as in Häffner .
(2023). The large language model we use is ConfliBERT (Y. Hu ., 2022).
ConfliBERT has the same architecture as BERT-base and is pretrained from
scratch using a large corpus in the politics and conflicts
domain.141414https://github.com/eventdata/ConfliBERT/. When we finetune
ConfliBERT using CrisisWatch texts, we set the learning rate to 2e-5 and set
the number of training epochs from 10 (see Listing 2).

⬇

1BASE_MODEL = "snowood1/ConfliBERT-scr-uncased" # pretrained model

2LEARNING_RATE = 2e-5

3MAX_LENGTH = 256

4BATCH_SIZE = 16

5EPOCHS = 10

6

7training_args = TrainingArguments(

8 learning_rate=LEARNING_RATE, # learning rate

9 per_device_train_batch_size=BATCH_SIZE, # training batch size

10 per_device_eval_batch_size=BATCH_SIZE, # test batch size

11 num_train_epochs=EPOCHS, # number of training epochs

12 seed=123

13)

14

15tokenizer = AutoTokenizer.from_pretrained(BASE_MODEL)

16

17def model_init():

18 model = AutoModelForSequenceClassification.from_pretrained(BASE_MODEL,
num_labels=1, ignore_mismatched_sizes=True) # number of labels

19 return model

20

21def preprocess_function(examples):

22 label = examples["fatalities_log"]

23 examples = tokenizer(examples["final_text"], truncation=True,
padding="max_length", max_length=256) # max sequence length

24 examples["label"] = float(label)

25 return examples

26

27class RegressionTrainer(Trainer):

28 def compute_loss(self, model, inputs, return_outputs=False):

29 labels = inputs.pop("labels")

30 outputs = model(**inputs)

31 logits = outputs[0][:, 0]

32 loss = torch.nn.functional.mse_loss(logits, labels)

33 return (loss, outputs) if return_outputs else loss

Listing 2: Specify the number of labels to be 1 for the regression task and
set the hyperparameters such as batch size and learning rate.

We report two groups of finetuning results: ConfliBERT Length-256, where we
set the max number of tokens to 256 and ConfliBERT Length-512, where we set
the maximum sequence length to 512 (Wang, 20232).151515We note that 3,386 out
of 23,868 samples (14%) contain more than 256 tokens. By setting the max
sequence length to 256, we are effectively truncating these long samples to
256 tokens, thus reducing the amount of information that we give to the model.
Setting the max sequence length to 512, which is the longest input length
possible for BERT models, helps alleviate this problem. In Table 2, we compare
the performance of finedtuned ConfliBERT models with that of two baseline
models. OCoDi-Random Forest is the dictionary-based model that leverages
random forest, where OCoDi stands for Objective Conflict Dictionary. OCoDi-
XGBoost is the dictionary-based model that leverages XGBoost (Häffner .,
2023).

Table 2: Finetuned ConfliBERT models outperform dictionary-based models by a large margin. Results in Columns 1 and 2 are from Table 2 in Häffner . (2023). By increasing the maximum sequence length to 512, we are able to further improve the performance of those finetuned models. Best results in bold. | Model
---
| OCoDi
---
Random Forest
(1)
| OCoDi
---
XGBoost
(2)
| ConfliBERT
---
Length-256
(3)
| ConfliBERT
---
Length-512
(4)
MSE | 1.59 | 1.60 | 0.99 | 0.82
R2 | 0.64 | 0.63 | 0.77 | 0.81

We observe that by finetuning ConfliBERT (Column 3), we are able to achieve a
much lower mean squared error and a substantially higher R2 than the
dictionary-based approaches. Further, by increasing the maximum sequence
length from 256 to 512 (Column 4), we observe that the mean squared error on
the test set decreases from 0.99 to 0.82 and that the R2 on the test set
increases from 0.77 to 0.81. Comparing OCoDi-XGBoost and ConfliBERT Max
Length, we observe that a finetuned ConfliBERT is able to achieve an MSE that
is 49% lower than OCoDi-XGBoost and an R2 that is 29% higher.

##  5 Conclusion

The pretrain-finetune paradigm has revolutionized the field of natural
language processing. In this tutorial, we have provided an intuitive and
thorough walk-through of the key concepts therein. We also introduced some of
the most common use cases that the pretrain-finetune paradigm supports. To
help facilitate the wider adoption of this new paradigm, we have included
easy-to-follow examples and made the related datasets and scripts publicly
available. Practitioners and scholars in quantitative psychology and
psychology more broadly should find our tutorial useful.

## References

  * Brown . (2020) gpt3Brown, T., Mann, B., Ryder, N., Subbiah, M., Kaplan, JD., Dhariwal, P.Amodei, D.  2020\.  Language Models are Few-Shot Learners Language models are few-shot learners.  34th Conference on Neural Information Processing Systems (NeurIPS 2020).
  * Devlin . (2019) bertDevlin, J., Chang, MW., Lee, K. Toutanova, K.  2019\.  BERT: Pre-training of Deep Bidirectional Transformers for Language Understanding BERT: Pre-training of Deep Bidirectional Transformers for Language Understanding.  Proceedings of NAACL-HLT4171-4186.
  * E. Hu . (2022) loraHu, E., Shen, Y., Wallis, P., Allen-Zhu, Z., Li, Y., Wang, S.Chen, W.  2022\.  LoRA: Low-Rank Adaptation of Large Language Models Lora: Low-rank adaptation of large language models.  Proceedings of the International Conference on Learning Representations.
  * Y. Hu . (2022) conflibertHu, Y., Hosseini, M., Parolin, ES., Osorio, J., Khan, L., Brandt, PT. D’Orazio, VJ.  2022\.  ConfliBERT: A Pre-trained Language Model for Political Conflict and Violence Conflibert: A pre-trained language model for political conflict and violence.  Proceedings of the 2022 Conference of the North American Chapter of the Association for Computational Linguistics.
  * Häffner . (2023) interpretableHäffner, S., Hofer, M., Nagl, M. Walterskirchen, J.  2023\.  Introducing an Interpretable Deep Learning Approach to Domain-Specific Dictionary Creation: A Use Case for Conflict Prediction Introducing an interpretable deep learning approach to domain-specific dictionary creation: A use case for conflict prediction.  Political Analysis1–19.  10.1017/pan.2023.7
  * Kjell . (2023) text_packageKjell, O., Giorgi, S. Schwartz, HA.  2023\.  The Text-Package: An R-Package for Analyzing and Visualizing Human Language Using Normal Language Processing and Transformers The text-package: An r-package for analyzing and visualizing human language using normal language processing and transformers.  Psychological Methods.
  * Lan . (2020) albertLan, Z., Chen, M., Goodman, S., Gimpel, K., Sharma, P. Soricut, R.  2020\.  ALBERT: A Lite BERT for Self-supervised Learning of Language Representations ALBERT: A Lite BERT for Self-supervised Learning of Language Representations.  ICLR.
  * Liu . (2019) robertaLiu, Y., Ott, M., Goyal, N., Du, J., Joshi, M., Chen, D.Stoyanov, V.  2019\.  RoBERTa: A Robustly Optimized BERT Pretraining Approach RoBERTa: A Robustly Optimized BERT Pretraining Approach.  arXiv:1907.11692.
  * Longpre . (2020) augLongpre, S., Wang, Y. DuBois, C.  2020\.  How Effective is Task-Agnostic Data Augmentation for Pretrained Transformers? How Effective is Task-Agnostic Data Augmentation for Pretrained Transformers?  Findings of the Association for Computational Linguistics: EMNLP 2020.
  * Mikolov . (2013) word2vecMikolov, T., Sutskever, I., Chen, K., Corrado, G. Dean, J.  2013\.  Distributed Representations of Words and Phrases and their Compositionality Distributed representations of words and phrases and their compositionality.  NIPS’13: Proceedings of the 26th International Conference on Neural Information Processing Systems.
  * Osnabrügge . (2021) cross_domainOsnabrügge, M., Ash, E. Morelli, M.  2021\.  Cross-Domain Topic Classification for Political Texts Cross-Domain Topic Classification for Political Texts.  Political Analysis.
  * Pargent . (2023) best_practices_supervised_mlPargent, F., Schoedel, R. Stachl, C.  2023\.  Best Practices in Supervised Machine Learning: A Tutorial for Psychologists Best practices in supervised machine learning: A tutorial for psychologists.  Advances in Methods and Practices in Psychological Science.
  * Vaswani . (2017) attentionVaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, AN.Polosukhin, I.  2017\.  Attention Is All You Need Attention Is All You Need.  31st Conference on Neural Information Processing Systems.
  * Wang (2019) predicted_acc2Wang, Y.  2019\.  Comparing Random Forest with Logistic Regression for Predicting Class-Imbalanced Civil War Onset Data: A Comment Comparing Random Forest with Logistic Regression for Predicting Class-Imbalanced Civil War Onset Data: A Comment.  Political Analysis211107-110.
  * Wang (20231) finetune_paWang, Y.  20231\.  On Finetuning Large Language Models On finetuning large language models.  Political Analysis.
  * Wang (20232) pretrained_topic_classificationWang, Y.  20232\.  Topic Classification for Political Texts with Pretrained Language Models Topic Classification for Political Texts with Pretrained Language Models.  Political Analysis.
  * Wang . (2022) bag_of_wordsWang, Y., Tian, J., Yazar, Y., Ones, DS. Landers, RN.  2022\.  Using Natural Language Processing and Machine Learning to Replace Human Content Coders Using natural language processing and machine learning to replace human content coders.  Psychological Methods.
  * Zhang . (2021) monitor_depressionZhang, Y., Lyu, H., Liu, Y., Zhang, X., Wang, Y. Luo, J.  2021\.  Monitoring Depression Trends on Twitter During the COVID-19 Pandemic: Observational Study Monitoring depression trends on twitter during the covid-19 pandemic: Observational study.  JMIR Infodemiology.