# Feature-Based Time-Series Analysis in R using the theft Package

Trent Henderson
The University of Sydney &Ben D. Fulcher
The University of Sydney<EMAIL_ADDRESS>
Trent Henderson, Ben Fulcher The theft Package for R Time series are measured
and analyzed across the sciences. One method of quantifying the structure of
time series is by calculating a set of summary statistics or ‘features’, and
then representing a time series in terms of its properties as a feature
vector. The resulting feature space is interpretable and informative, and
enables conventional statistical learning approaches, including clustering,
regression, and classification, to be applied to time-series datasets. Many
open-source software packages for computing sets of time-series features exist
across multiple programming languages, including catch22 (22 features: Matlab,
R, Python, Julia), feasts (42 features: R), tsfeatures (63 features: R), Kats
(40 features: Python), tsfresh (779 features: Python), and TSFEL (390
features: Python). However, there are several issues: (i) a singular access
point to these packages is not currently available; (ii) to access all feature
sets, users must be fluent in multiple languages; and (iii) these feature-
extraction packages lack extensive accompanying methodological pipelines for
performing feature-based time-series analysis, such as applications to time-
series classification. Here we introduce a solution to these issues in the
form of a statistical software package for R called theft: Tools for Handling
Extraction of Features from Time series. theft is a unified and extendable
framework for computing features from the six open-source time-series feature
sets listed above. It also includes a suite of functions for processing and
interpreting the performance of extracted features, including extensive data-
visualization templates, low-dimensional projections, and time-series
classification operations. With an increasing volume and complexity of large
time-series datasets in the sciences and industry, theft provides a
standardized framework for comprehensively quantifying and interpreting
informative structure in time series. time-series analysis, time-series
features, R, machine learning time-series analysis, time-series features, R,
machine learning Trent Henderson
Dynamics and Neural Systems Group
School of Physics
The University of Sydney
Camperdown NSW 2006, Australia
E-mail:

##  1 Introduction

Taking repeated measurements of some quantity through time, forming a time
series, is common across the sciences and industry. The types of time series
commonly analyzed are diverse, ranging from signals from an
electroencephalogram (West _et al._ , 1999), CO2 concentration in the
atmosphere (Kodra _et al._ , 2011), light-curves from distant stars (Barbara
_et al._ , 2022), and the number of clicks on a webpage (Kao _et al._ , 2021).
We can ask many different questions about such data, for example: (i) “can we
distinguish the dynamics of brain disorders from neurotypical brain
function?”; (ii) “can we classify different geospatial regions based on their
temporal CO2 concentration”; or (iii) “can we classify new stars based on
their light curves?”. One approach to answering such questions is to capture
properties of each time series and use that information to train a
classification algorithm. This can be achieved by extracting from each time
series a set of interpretable summary statistics or ‘features’. Using this
procedure, a collection of univariate time series can be represented as a time
series $\times$ feature matrix which can be used as the basis for a range of
conventional statistical learning procedures (Fulcher _et al._ , 2013;
Fulcher, 2018).

The range of time-series analysis methods that can be used to define time-
series features is vast, including properties of the distribution,
autocorrelation function, stationarity, entropy, methods from the physics
nonlinear time-series analysis literature (Fulcher _et al._ , 2013). Because
features are direct outputs of a mathematical operation, and are often tightly
linked to underlying theory (e.g., Fourier analysis or information theory),
they can yield interpretable understanding of patterns in time series and the
processes that produce them—information that can guide further investigation.
The first work to organize these methods from across the interdisciplinary
literature encoded thousands of diverse time-series analysis methods as
features and compared their behavior on a wide range of time series (Fulcher
_et al._ , 2013). The resulting interdisciplinary library of thousands of
time-series features has enabled new ways of doing time-series analysis,
including the ability to discover high-performing methods for a given problem
in a systematic, data-driven way through large-scale comparison (overcoming
the subjective and time-consuming task of selecting methods manually) (Fulcher
and Jones, 2014). This approach has been termed ‘highly comparative time-
series analysis’, and has been implemented in the Matlab software hctsa, which
computes $>7700$ time-series features (Fulcher and Jones, 2017). The approach
of automated discovery provided by hctsa has been applied successfully to many
scientific problems, such as classifying zebra finch motifs across different
social contexts (Paul _et al._ , 2021), classifying cord pH from fetal heart-
rate dynamics (Fulcher _et al._ , 2012), and classifying changes in cortical
dynamics from manipulating the firing of excitatory and inhibitory neurons
(Markicevic _et al._ , 2020). While hctsa is comprehensive in its coverage of
time-series analysis methods, calculating all of its features on a given
dataset is computationally expensive and it requires access to the proprietary
Matlab software, limiting its broader use.

The past decade has seen the development of multiple software libraries that
implement different sets of time-series features across a range of open-source
programming languages. Here, we focus on the following six libraries:

  * •

catch22 (C, Matlab, R, Python, Julia) computes a representative subset of 22
features from hctsa (Lubba _et al._ , 2019). The $>$7700 features in hctsa
were applied to 93 time-series classification tasks to retain the smallest
number of features that maintained high performance on these tasks while also
being minimally redundant with each other, yielding the catch22 set. catch22
was coded in C for computational efficiency, with wrappers for Matlab, and
packages for: R, as Rcatch22 (Henderson, 2021); Julia, as Catch22.jl (Harris,
2021); and Python, as pycatch22. The construction of the catch22 feature set
focused on dynamical properties, but users can also include the mean and
standard deviation of time series in addition to the regular 22 features
(yielding the 24-feature set termed ‘catch24’).

  * •

tsfeatures (R) is the most prominent package for computing time-series
features in R (Hyndman _et al._ , 2020). The 63 features in tsfeatures include
techniques commonly used by econometricians and forecasters, such as crossing
points, seasonal and trend decomposition using Loess (STL; Cleveland _et al._
(1990)), autoregressive conditional heteroscedasticity (ARCH) models, unit-
root tests, and sliding windows. tsfeatures also includes a small subset of
features from hctsa that were previously used to organize tens of thousands of
time series in the CompEngine time-series database (Fulcher _et al._ , 2020).

  * •

feasts (R) shares a subset of the same features as tsfeatures, computing a
total of 42 features. However, the scope of feasts as a software package is
larger: it is a vehicle to incorporate time-series features into the software
ecosystem known as the tidyverts111https://tidyverts.org—a collection of
packages for time series that follow tidy data principles (Wickham, 2014).
This ensures alignment with the broader and popular tidyverse collection of
packages for data wrangling, summarization, and statistical graphics (Wickham
_et al._ , 2019). feasts also includes functions for producing graphics, but
these are largely focused on exploring quantities of interest in econometrics,
such as autocorrelation, seasonality, and STL decomposition.

  * •

tsfresh (Python) includes 779 features that measure properties of the
autocorrelation function, entropy, quantiles, fast Fourier transforms, and
distributional characteristics (Christ _et al._ , 2017). tsfresh also includes
a built-in feature filtering procedure, FeatuRe Extraction based on Scalable
Hypothesis tests (FRESH), that uses a hypothesis-testing process to control
the percentage of irrelevant extracted features (Christ _et al._ , 2018).
tsfresh has been used widely to solve time-series problems, such as anomaly
detection in Internet-of-Things streaming data (Yang _et al._ , 2021) and
sensor-fault classification (Liu _et al._ , 2020).

  * •

TSFEL (Python) contains 390 features that measure properties associated with
distributional characteristics, the autocorrelation function, fast Fourier
transforms, spectral quantities, and wavelets (Barandas _et al._ , 2020).
TSFEL was initially designed to support feature extraction of inertial
data—such as data produced by human wearables—for the purpose of activity
detection and rehabilitation.

  * •

Kats (Python), developed by Facebook Research, contains a broad range of time-
series functionality, including operations for forecasting, outlier and
property detection, and feature calculation (Facebook Infrastructure Data
Science, 2021). The feature-calculation module of Kats is called TSFeatures
and includes 40 features (30 of which are based on R’s tsfeatures package).
Kats includes features associated with crossing points, STL decomposition,
sliding windows, autocorrelation and partial autocorrelation, and Holt–Winters
methods for detecting linear trends.

The six sets vary over several orders of magnitude in their computation time,
and exhibit large differences in both within-set feature redundancy—how
correlated features are within a given set—and between-set feature
redundancy—how correlated, on average, features are between different pairwise
comparisons of sets (Henderson and Fulcher, 2021). While each set contains a
range of features that could be used to tackle time-series analysis problems,
there are currently no guidelines for selecting an appropriate feature set for
a given problem, nor methods for combining the different strengths of all
sets. Performance on a given time-series analysis task depends on the choice
of the features that are used to represent the time series, highlighting the
importance of being able to easily compute many different features from across
different feature sets. Furthermore, following feature extraction, there is no
set of visualization and analysis templates for common feature-based problem
classes, such as feature-based time-series classification (like the tools
provided in hctsa Fulcher and Jones (2017)). Here we present a solution for
these challenges in the form of an open-source package for R called theft:
Tools for Handling Extraction of Features from Time series.

##  2 The theft package for R

theft unifies the six free and open-source feature sets described in Section
1, thus overcoming barriers in using diverse feature sets developed in
different software environments and using different syntax. theft also
provides an extensive analytical pipeline and statistical data visualization
templates similar to those found in hctsa for understanding feature behavior
and performance. Such pipelines and templates do not currently exist in the
free and open-source setting, making theft an invaluable tool for both
computing and understanding features. While there is some software support for
computing features in a consistent setting (such as in tsflex Van Der Donckt
_et al._ (2022), which also provides sliding window extraction capability),
such software is limited to specifying the functional form of individual time-
series features rather than automatically accessing every feature contained in
different sets.

Figure 1:  theft implements a workflow for extracting features from univariate
time series and processing and analyzing the results. First, a time-series
dataset (A) is converted into a tidy (‘long’) data frame (B) with variable
names for unique identifiers, time-point indices, values, and group labels
(e.g., in the case of classification problems). If one of the feature sets
selected is a Python library, the user can point R to the Python version
containing the installed software (C). One or more feature sets are then
computed on the dataset (D). A range of statistical analysis and data
visualization functionality is also implemented, including: (E) feature
quality assessment (e.g., understanding the proportion of non-NA values by
feature); (F) normalized time series $\times$ feature matrix visualization;
(G) low-dimensional projections of the feature space; and (H) normalized
feature $\times$ feature correlation matrix visualization. Functionality is
also provided for time-series classification (a common application of feature-
based time-series analysis), including: (I) understanding the most
discriminative individual features; and (J) fitting and evaluating classifiers
with more than one feature as input.

The functionality provided by theft is summarized in Fig. 1 and broadly
follows the feature-based time-series analysis workflow of hctsa (Fulcher and
Jones, 2017). The workflow begins with a time-series dataset (Fig. 1A) that is
converted to a tidy format (Fig. 1B). If any of the Python feature sets are to
be used, the Python environment containing the installed software is passed to
theft using init_theft (Fig. 1C). Time-series features are then extracted
(Fig. 1D). The user can pass the extracted features into a range of
statistical and visualization functions to derive interpretable understanding
of the informative patterns in their dataset (Figs 1E–J). Importantly, theft
uses R’s S3 object-oriented programming system, meaning classes and their
methods are defined to ensure usability. Classes are defined for feature
calculation objects (D; feature_calculations) and low dimensional projection
objects (G; low_dimension). Almost every theft function takes at least an
object of one of these classes as input. This means that generic methods, such
as plot can be called on the object, where additional arguments to the generic
function such as type declares the type of graphic that is produced. These
methods are discussed in detail in the following sections.

In this paper, we demonstrate how theft can be used to tackle a time-series
classification problem, using the Bonn University electroencephalogram (EEG)
dataset as a case study (Andrzejak _et al._ , 2001). The dataset contains 500
time series, with 100 time series each from five labeled classes: (i) awake
with eyes open (labeled ‘eyesOpen’); (ii) awake with eyes closed
(‘eyesClosed’); (iii) epileptogenic zone (‘epileptogenic’); (iv) hippocampal
formation of the opposite hemisphere of the brain (‘hippocampus’); and (v)
seizure activity (‘seizure’). Note that classes (i) and (ii) are from healthy
volunteers, while classes (iii), (iv), and (v) are from a presurgical
diagnosis archive. This dataset was chosen as a demonstrative example because
it has been widely studied as a time-series classification problem, and prior
studies have focused on properties of the dynamics that accurately distinguish
the classes—which is well-suited to the feature-based approach. For example,
an analysis using hctsa revealed that seizure recordings are characterized
most notably by higher variance, as well as lower entropy, lower long-range
scaling exponents, and many other differences (Fulcher _et al._ , 2013).
Further, it was also found that 172 individual features within hctsa could
distinguish between healthy EEGs and seizures using a 10-fold cross-validation
linear classifier with $>95\%$ accuracy, with eight of these features
achieving $>98.75\%$ — exceeding previous results which used operations from
the discrete wavelet transform as inputs to a support vector machine
classifier (Subasi and Ismail Gursoy, 2010).

###  2.1 Extracting features

In feature-based time-series analysis, each univariate time series in a
dataset is represented as a feature vector, such that the dataset can be
represented as a time series $\times$ feature data matrix. Any single feature
set, or combination of multiple feature sets, can be computed for a given
time-series dataset with the theft function calculate_features. An example
call which extracts features from all six sets is shown in Listing 1. Note
that throughout this paper we will present the code to execute each piece of
analysis, but will largely omit the extensive amount of optional arguments for
clarity. The vignette and function documentation that comes with every
download of theft contains detailed information regarding these arguments
(Henderson and Bryant, 2022). The call in Listing 1 takes a tidy data frame
that contains the time-series data (tmp in this example), takes the relevant
user-specified columns, and computes the time series $\times$ feature matrix
for the specified feature set(s). The output of this function is an S3 object
of class feature_calculations, which in this example is stored in the R
environment as feature_matrix. Within this object is a data frame which
contains five columns if the dataset is labeled (as in time-series
classification), and four otherwise: id (unique identifier for each time
series), names (feature name), values (feature value), method (feature set),
and group (class label, if applicable). This output structure ensures that,
regardless of the feature set selected, the resulting object is always of the
same format and can be used with the rest of theft’s functions without manual
data reshaping. For the Bonn EEG dataset, calculating features for all sets in
theft takes several hours.

⬇

1 all_features <- calculate_features(

2 data = tmp,

3 id_var = "id",

4 time_var = "timepoint",

5 values_var = "values",

6 group_var = "group",

7 feature_set = c("catch22", "feasts", "tsfeatures",

8 "tsfresh", "TSFEL", "Kats")

Listing 1:

###  2.2 Assessing feature extraction quality

Not all features return real-valued outputs for all time series, meaning non-
numeric values, such as NaN or Inf/-Inf, or even missing outputs or errors,
can be returned as a result of feature extraction. For example, a feature
which computes the variance across multiple 200-sample time-series windows
cannot be computed for time series shorter than 200 samples (as there are not
enough samples to form even a single window). For effective quality control,
it is important to visualize numeric and non-numeric outputs following feature
extraction, which is implemented in theft by calling plot on a
feature_calculations object and specifying a plot type of "quality", as shown
in Listing 2. This graphic plots the proportion of values in each that are
numeric, NaN/NA, or Inf/-Inf for each feature as bars.

⬇

1 plot(all_features, type = "quality")

Listing 2:

###  2.3 Normalizing features

Different features vary over very different ranges; e.g., features that
estimate $p$-values from a hypothesis test vary over the unit interval,
whereas a feature that computes the length of a time series can take (often
large) positive integers. These differences in scale can complicate the
visualization of feature behavior and the construction of statistical learning
algorithms involving diverse features. To overcome these limitations, a common
pre-processing step involves scaling all features. theft includes four such
methods for converting a set of raw feature values, $\mathbf{x}$, to a
normalized version, $\mathbf{z}$:

  1. 1.

$z$-score: $z_{i}=\frac{x_{i}-\mu}{\sigma}$,

  2. 2.

linear scaling to unit interval:
$z_{i}=\frac{x_{i}-\min(\mathbf{x})}{\max(\mathbf{x})-\min(\mathbf{x})}$,

  3. 3.

sigmoid: $z_{i}=\left[1+\exp(-\frac{x_{i}-\mu}{\sigma})\right]^{-1}$,

  4. 4.

and outlier-robust sigmoid:
$z_{i}=\left[1+\exp\left(-\frac{x_{i}-\mathrm{median}(\mathbf{x})}{\mathrm{IQR}(\mathbf{x})/{1.35}}\right)\right]^{-1}$,

where $\mu$ is the mean, $\sigma$ is the standard deviation, and
$\mathrm{IQR}(\mathbf{x})$ is the interquartile range of $\mathbf{x}$. All
four transformations end with a linear rescaling to the unit interval. The
outlier-robust sigmoid transformation, introduced in Fulcher _et al._ (2013),
can be helpful in normalizing feature-value distributions with large outliers.
Normalization is an option in each of the core analysis and visualization
functions within theft, but users can also perform normalization outside of
these functions on individual feature vectors or feature_calculations objects
(i.e., the entire feature matrix) using the normalise function which
automatically detects the input object class. An example for
feature_calculations objects is shown in Listing 3.

⬇

1 normed <- normalise(all_features, method = "RobustSigmoid")

Listing 3:

###  2.4 Visualizing the feature matrix

A hallmark of large-scale feature extraction is the ability to visualize the
intricate patterns of how different scientific algorithms behave across a
time-series dataset. This can be achieved in theft using by specifying type =
"matrix" when calling plot on a feature_calculations object to produce a
heatmap of the time series (rows) $\times$ feature matrix (columns) which
organizes the rows and columns to help reveal interesting patterns in the
data. The plot of the combination of all six open feature sets for the Bonn
EEG dataset is shown in Fig. 2, with the code displayed in Listing 4. We can
see some informative structure in this graphic, including many groups of
features with similar behavior on this dataset (i.e., columns with similar
patterns), indicating substantial redundancy across the joint set of features
(Henderson and Fulcher, 2021). The top block of 100 rows, which visually have
the most distinctive properties, were found to correspond to time series from
the “seizure” class, indicating the ability of this large combination of time-
series features to meaningfully structure the dataset.

Figure 2:  A time series $\times$ feature matrix heatmap produced by
generating a matrix plot on the feature_calculations object. Extracted feature
vectors for each time series (500) in the Bonn EEG dataset using all six
feature sets in theft (1253 features in total, after filtering out 63 features
with NaN values) are represented as a heatmap. Similar features (columns) and
time series (rows) are positioned close to each other using (average)
hierarchical clustering. Each tile is a normalized value for a given time
series and feature. This plot is generated by the code in Listing 4.

In matrix plots in theft, hierarchical clustering is used to reorder rows and
columns so that time series (rows) with similar properties are placed close to
each other and features (columns) with similar behavior across the dataset are
placed close to each other—where similarity in behavior is quantified using
Euclidean distance in both cases (Day and Edelsbrunner, 1984).

Default settings within plot enable users to easily generate outputs in a
single line of code, but more advanced users may seek to tweak the optional
arguments. For example, different linkage algorithms for hierarchical
clustering can be controlled supplied to clust_method, which uses average
linkage as a default, and different rescaling methods can be supplied to the
method argument, which defaults to "z-score".

⬇

1 plot(all_features, type = "matrix")

Listing 4:

###  2.5 Projecting low-dimensional feature-spaces

Low-dimensional projections are a useful tool for visualizing the structure of
high-dimensional datasets in low-dimensional spaces. Here we are interested in
representing a time-series dataset in a two-dimensional projection of the
high-dimensional feature space, which can reveal structure in the dataset,
including how different labeled classes are organized. For linear
dimensionality reduction techniques—such as principal components analysis
(PCA) (Jolliffe, 2002)—the results can be visualised in two dimensions as a
scatterplot, where the principal component (PC) that explains the most
variance in the data is positioned on the horizontal axis and the second PC on
the vertical axis, and each time series is represented as a point (colored by
its group label in the case of a labeled dataset). When the structure of a
dataset in the low-dimensional feature space matches known aspects of the
dataset (such as class labels), it suggests that the combination of diverse
time-series features can capture relevant dynamical properties that differ
between the classes. It can also reveal new types of structure in the dataset,
like clear sub-clusters within a labeled class, that can guide new
understanding of the dataset. Low-dimensional projections of time-series
features have been shown to meaningfully structure time-series
datasets—revealing sex and day/night differences in Drosophila (Fulcher and
Jones, 2017), distinguishing types of stars based on their light curves
(Barbara _et al._ , 2022), and categorizing sleep epochs (Decat _et al._ ,
2022).

In theft, both a linear dimensionality reduction method—PCA—and a nonlinear
dimensionality reduction method—$t$-distributed stochastic neighbour embedding
($t$-SNE) (van der Maaten and Hinton, 2008)—are included. While many
dimensionality-reduction algorithms exist (Sorzano _et al._ , 2014), here we
selected just these two to minimize package dependencies. Time-series datasets
can be projected in low-dimensional feature spaces in theft using the
reduce_dims function demonstrated in Listing 5. For $t$-SNE, users can control
the perplexity hyperparameter using the perplexity argument. The reduce_dims
function returns an S3 object of class low_dimension, which can then be passed
to plot which will automatically detect the class. The plot method for this
object contains only one other argument (show_covariance) which is a Boolean
specifying whether to draw covariance ellipses for each group in the
scatterplot (if a grouping variable is detected).

⬇

1 low_dim_calc <- reduce_dims(all_features, method = "MinMax", low_dim_method
= "t-SNE", perplexity = 15)

Listing 5:

The low-dimensional projection plot for the Bonn EEG dataset (using $t$-SNE
and all non-NaN features across the six feature sets included in theft) is
shown in Fig. 3 with perplexity 10, as produced by the code in Listing 5. The
low-dimensional projection (formed from $>1200$ features in theft)
meaningfully structures the labeled classes of the dataset. Specifically, two
of the presurgical diagnosis classes—“epileptogenic” (epileptogenic zone) and
“hippocampus” (hippocampal formation of the opposite hemisphere of the
brain)—appear to exhibit considerable overlap in the projected space, while
the two healthy volunteer classes “eyesOpen” (awake state with eyes open) and
“eyesClosed” (awake state with eyes closed) occupy space further away from the
other classes but closer to each other. The “seizure” class occupies a space
largely separate from the other four classes in the projection, consistent
with its distinctive dynamics (Fulcher _et al._ , 2013).

Figure 3:  Low-dimensional projection of the Bonn EEG dataset using theft.
Using $t$-SNE with perplexity 15, the high-dimensional feature space of over
$1200$ features is projected into two dimensions. Each point represents a time
series which is colored according to its class label. Time series that are
located close in this space have similar properties, as measured by the six
feature sets in theft. This plot is generated by calling plot(low_dim_calc,
show_covariance = FALSE) after the code in Listing 5.

###  2.6 Constructing classifiers with multiple features

Combinations of complementary, discriminative features can often be used to
construct accurate time-series classifiers (Fulcher and Jones, 2014). Drawing
on computed time-series features (that may derive from one or more existing
feature sets), theft can fit and evaluate classifiers using the
fit_multi_feature_classifier function. This allows users to evaluate the
relative performance of each feature set, of the combination of all sets, or
any other combination of features. Providing easy access to a range of
classification algorithms and accompanying inferential tools (such as null
permutation testing to obtain $p$-values) through fit_multi_feature_classifier
allows users to compare sets of features to better understand the most
accurate feature sets for a given time-series classification problem. The code
presented in Listing 6 provides an example usage for the Bonn EEG dataset with
a linear support vector machine (SVM) classifier.

⬇

1 mf_results <- fit_multi_feature_classifier(

2 data = all_features,

3 by_set = TRUE,

4 test_method = "svmLinear",

5 use_k_fold = TRUE,

6 num_folds = 10)

Listing 6:

The fit_multi_feature_classifier function returns a list object. If by_set is
TRUE, an plot object called FeatureSetResultsPlot is created and returned in
the list which contains a bar plot of classification accuracy for each feature
set (if by_set is set to FALSE, all features will be used as predictors in the
chosen classification model, ignoring the set they originate from and not
returning a plot).

The modeling components of fit_multi_feature_classifier are executed through a
wrapper for the machine-learning package caret (Kuhn, 2020). This means the
extensive list of classification models available in caret can be accessed
from theft by specifying the method name in the test_method argument. Prior to
fitting a model, fit_multi_feature_classifier performs two operations: (i)
filtering out features that are constants or contain NaN/NA or Inf/-Inf
values; and (ii) re-coding of class labels into syntactically valid names for
model objects in R. The resulting data is then passed into a caret train
operation, where if use_k_fold is set to TRUE, $k$-fold cross-validation is
performed, with the number of folds set by the num_folds argument. All
operations produced by fit_multi_feature_classifier use centering and scaling
preprocessing provided by caret and exclude time-series features with near-
zero variance after executing this procedure. The performance metric is
specified by use_balanced_accuracy where if FALSE, mean classification
accuracy is used and if TRUE, balanced mean classification accuracy is used.
Balanced classification accuracy is a useful metric for problems where class
imbalances can artificially inflate the accuracy metric (i.e., a classifier
might assign the majority or all of the predicted class labels to the class
with the most representation in the data).

The summary FeatureSetResultsPlot graphic produced by Listing 6 for the Bonn
EEG dataset (in which all time series have been $z$-scored, to focus on
differences in dynamical properties between the classes) is shown in Fig. 4.
On this problem, we find that TSFEL (with 378 features after filtering) has
the highest mean classification accuracy ($70.8\%$, SD = $3.6\%$ over 10
folds) and feasts has the lowest mean accuracy ($60.4\%$, SD = $4.2\%$). The
combination of all $>1260$ features from across every set demonstrates
accuracy consistent with the top performer (TSFEL), with slightly lower mean
performance (potentially due to over-fitting in a higher-dimensional feature
space).

Figure 4:  Comparison of mean classification accuracy between feature sets in
theft for the five-class Bonn EEG classification task. Classification accuracy
using a linear SVM is presented for each of the six feature sets in theft as
well as the combination of all their features (‘All features’). The number of
features retained for analysis after filtering is displayed in parentheses
after the feature set name on the horizontal axis which has been sorted from
highest to lowest mean accuracy. Mean classification accuracy across the same
10 cross-validation folds is displayed as colored points for each set with
$\pm 1$SD error bars. This plot is generated by the code in Listing 6.

####  2.6.1 Permutation testing

In applications involving small datasets, or when small effects are expected,
it is useful to quantify how different the calculated classification
performance is from a null setting in which data are classified randomly. One
method for inferring test statistics is to use permutation testing—a procedure
that samples a null process many times to form a distribution against which a
value of importance (i.e., the classification accuracy result from a model)
can be compared to estimate a $p$-value (Ojala and Garriga, 2009). In theft,
permutation testing is implemented for evaluating classification performance
in fit_multi_feature_classifier through the use_empirical_null argument. When
set to TRUE, an object is returned in the list called TestStatistics which
contains a data frame of classification accuracy results and associated
$p$-values. Regardless of the options specified to by_set and
use_empirical_null, an object called RawClassificationResults is always
returned in the overall list object created by fit_multi_feature_classifier,
which contains a data frame of classification accuracy outputs from each model
that is trained and evaluated.

A more detailed version of Listing 6 with the optional arguments specified to
execute the null testing procedure is displayed in Listing 7.

⬇

1 mf_results <- fit_multi_feature_classifier(

2 data = all_features,

3 by_set = TRUE,

4 test_method = "svmLinear",

5 use_empirical_null = TRUE,

6 null_testing_method = "ModelFreeShuffles",

7 p_value_method = "gaussian",

8 num_permutations = 10000)

Listing 7:

Two methods for permutation testing are implemented: (i) model-free random
shuffles
("ModelFreeShuffles")—from the vector of class labels in the dataset, randomly
permutes num_permutations random shuffles (permutations) and calculates either
mean classification accuracy or balanced accuracy for each shuffled vector
against the original vector of class labels; and (ii) null model fits
("NullModelFits")—fits num_permutations models with the same classifier, same
input data, and same $k$-fold cross-validation procedure as the main
classification model, but trains the null models on randomly shuffled class
labels as the target variable.

The model-free procedure (which assumes a null model that produces randomized
outputs labels) is very fast, and provides a good approximation for the null
distribution obtained using randomized input labels when cross-validation is
used (i.e., when all evaluation is performed on unseen data). From the
resulting null distribution, a $p$-value can be assigned to the main model’s
performance statistic via one of two methods: (i) "gaussian" and (ii)
"empirical". The former takes the mean and standard deviation of the null
classification accuracy values, and evaluates the probability of the
classification accuracy result of the model fit on the correct class labels
against a Gaussian distribution parameterized by this null mean and standard
deviation. The latter simply calculates the proportion of null classification
accuracy values that are equal to or greater than the classification accuracy
result of the model fit on the correct class labels. If the null distribution
is believed to be approximately Gaussian, then an estimate of potentially
small $p$-values can be obtained using the "gaussian" setting with less
samples rather than expending a large computing time to resolve the potential
for small $p$-values through permutation testing and frequencies using the
"empirical" setting.

For the full Bonn EEG dataset, with 100 time series per class and strong
differences between signals, classification accuracies are far higher than
chance level (20%), and we obtain extremely small $p$-values, confirming the
low probability of obtaining such high accuracies by chance. To more clearly
demonstrate the permutation testing functionality, we analyzed a smaller
random subsample of the dataset: 14 $z$-scored time series each from the
‘hippocampus’ and ‘epileptogenic’ classes, and fit binary classification
models using the code presented in Listing 7. The mean classification accuracy
of feasts ($64\%$, $p=0.07$), catch22 ($68\%$, $p=0.07$), and TSFEL ($69\%$,
$p=0.07$) were not statistically significantly different (at $p<0.05$ with all
$p$-values adjusted using the Holm-Bonferroni method for multiple comparisons)
from the null distribution, despite their mean accuracy values being above the
chance probability of $50\%$. The remaining sets were statistically
significant: Kats ($74\%$, $p=0.02$), tsfresh ($79\%$, $p=0.007$), and
tsfeatures ($84\%$, $p=0.001$).

###  2.7 Finding and understanding informative individual features

Fitting models which use multiple features as inputs is often useful for
predicting class labels. However, users are also typically interested in
understanding patterns in their dataset, such as interpreting the types of
time-series analysis methods that best separate different classes, and the
relationships between these top-performing features. This can be achieved
using mass univariate statistical testing of individual features, quantifying
their importance either with conventional statistical tests (e.g., $t$-tests,
Wilcoxon Rank Sum Tests, and Signed Rank Tests), or with one-dimensional
classification algorithms (e.g., linear SVM, random forest classifiers). theft
implements the ability to identify top-performing features in the
compute_top_features function, with an example usage for the Bonn EEG dataset
(using features from all six packages) shown in Listing 8.

⬇

1 top_feature_results <- compute_top_features(

2 data = all_features,

3 num_features = 40,

4 test_method = "svmLinear",

5 use_k_fold = TRUE,

6 num_folds = 10,

7 use_empirical_null = TRUE)

Listing 8:

For two-class problems, users can access the traditional statistical tests by
specifying either test_method = "t-test", test_method = "wilcox", or
test_method =
"BinomialLogistic" to fit the desired statistical test instead of a caret
classification model. compute_top_features allows users to fit the same set of
caret classification models available in fit_multi_feature_classifier in the
one-dimensional space (i.e., the input to the algorithm is values on a single
time-series feature), which can be used for two-class or multi-class problems
(where traditional two-sample statistical tests cannot be used).

Regardless of the test_method used, compute_top_features always returns a list
object with three elements: (i) ResultsTable (a data frame of either
statistical test results and model information or classification results and
resulting null test statistics based on whether a caret model was used or
not); (ii) FeatureFeatureCorrelationPlot (pairwise feature $\times$ feature
correlation plot with the correlation method in the cor_method argument); and
(iii) ViolinPlots (violin plots for each feature that show each time series as
a point, colored and arranged in columns by class, with class-level
probability density lines around them).

The ResultsTable summary of classification accuracy values and $p$-values for
the top 40 features on the Bonn EEG classification task (determined by
$p$-value of mean classification accuracy against the empirical null) is
displayed in Table 1. There is representation from all six feature sets, with
the majority coming from the two largest sets: tsfresh and TSFEL. Inspecting
the feature names within the table can provide insight into the types of
features relevant for the classification problem. These features inlcude
properties associated with wavelets (e.g., 0_wavelet_energy_7, which measures
the area under the squared magnitude of the continuous wavelet transform at
scale 7), variance (e.g., 0_standard_deviation, which is the standard
deviation of the signal), and autocorrelation (e.g.,
values_autocorrelation_lag_7, which measures the value of the autocorrelation
function at lag 7). In addition, it can also be seen that two of strongest
performing features ($53.8\%$ individual classification accuracy) are in fact
the same feature implemented in two different sets
(SC_FluctAnal_2_rsrangefit_50_1_logi_prop_r1 from catch22 and
fluctanal_prop_r1 from tsfeatures). However, interpreting this table is
challenging as the relationships between the features are unknown—are all the
40 features behaving differently, or are they all highly correlated to each
other and essentially proxy metrics for the same underlying time-series
property? We can better understand these relationships by visualizing the
pairwise feature $\times$ feature correlation matrix.

Feature Feature Set Classification Accuracy $p$-value
SC_FluctAnal_2_rsrangefit_50_1_logi_prop_r1 catch22 53.8% $p$ $<$ .001
fluctanal_prop_r1 tsfeatures 53.8% $p$ $<$ .001 0_wavelet_energy_7 TSFEL 53.8%
$p$ $<$ .001 0_wavelet_standard_deviation_7 TSFEL 53.8% $p$ $<$ .001
0_wavelet_energy_5 TSFEL 52.0% $p$ $<$ .001 0_wavelet_standard_deviation_5
TSFEL 52.0% $p$ $<$ .001 0_wavelet_energy_4 TSFEL 51.4% $p$ $<$ .001
0_wavelet_standard_deviation_4 TSFEL 51.4% $p$ $<$ .001 0_wavelet_energy_8
TSFEL 51.0% $p$ $<$ .001 0_wavelet_standard_deviation_8 TSFEL 51.0% $p$ $<$
.001 values_autocorrelation_lag_6 tsfresh 50.6% $p$ $<$ .001
0_wavelet_energy_6 TSFEL 50.6% $p$ $<$ .001 0_wavelet_standard_deviation_6
TSFEL 50.6% $p$ $<$ .001 values_autocorrelation_lag_7 tsfresh 50.4% $p$ $<$
.001 seas_acf1 Kats 50.4% $p$ $<$ .001
values_agg_linear_trend_attr_stderr_chunk_len_10_f_agg_mean tsfresh 50.0% $p$
$<$ .001 hurst Kats 50.0% $p$ $<$ .001 0_wavelet_variance_6 TSFEL 49.8% $p$
$<$ .001 0_wavelet_variance_7 TSFEL 49.6% $p$ $<$ .001 0_spectral_distance
TSFEL 49.4% $p$ $<$ .001 values_autocorrelation_lag_5 tsfresh 48.8% $p$ $<$
.001 0_root_mean_square TSFEL 48.8% $p$ $<$ .001 0_area_under_the_curve TSFEL
48.8% $p$ $<$ .001 0_peak_to_peak_distance TSFEL 48.6% $p$ $<$ .001
values_agg_linear_trend_attr_stderr_chunk_len_5_f_agg_max tsfresh 48.4% $p$
$<$ .001 values_maximum tsfresh 48.4% $p$ $<$ .001
values_agg_linear_trend_attr_stderr_chunk_len_5_f_agg_mean tsfresh 48.4% $p$
$<$ .001 values_agg_linear_trend_attr_stderr_chunk_len_10_f_agg_max tsfresh
48.4% $p$ $<$ .001 values_agg_linear_trend_attr_stderr_chunk_len_50_f_agg_var
tsfresh 48.4% $p$ $<$ .001 0_wavelet_variance_8 TSFEL 48.4% $p$ $<$ .001 0_max
TSFEL 48.4% $p$ $<$ .001 0_mean_absolute_deviation TSFEL 48.2% $p$ $<$ .001
shift_level_max feasts 48.0% $p$ $<$ .001 var_tiled_mean feasts 47.8% $p$ $<$
.001 stability tsfeatures 47.8% $p$ $<$ .001
values_agg_linear_trend_attr_intercept_chunk_len_50_f_agg_min tsfresh 47.8%
$p$ $<$ .001 values_agg_linear_trend_attr_stderr_chunk_len_10_f_agg_min
tsfresh 47.8% $p$ $<$ .001 values_standard_deviation tsfresh 47.6% $p$ $<$
.001 values_linear_trend_attr_stderr tsfresh 47.6% $p$ $<$ .001
0_standard_deviation TSFEL 47.6% $p$ $<$ .001

Table 1:  Comparison of classification accuracy and $p$-values between the top
40 individual features in theft for the Bonn EEG dataset. $p$-values were
calculated by comparing each individual feature’s model classification
accuracy against a Gaussian null distribution parameterized by the mean and SD
of model-free shuffled samples. These values were generated by the code in
Listing 8.

The plot of the pairwise absolute correlation coefficients between the top 40
features (returned as FeatureFeatureCorrelationPlot in Listing 8) is displayed
in Fig. 5. The plot reveals two main groups of highly correlated
($|\rho|\gtrapprox 0.8$) features: in the bottom left and upper right of the
plot. The large cluster in the bottom left (containing features from tsfresh,
TSFEL, and feasts) contains features sensitive to signal variance. While a
number of features in this cluster have names associated with wavelets (e.g.,
TSFEL_0_wavelet_variance_7), the plot reveals that on this dataset, these
features exhibit similar behavior as features measuring signal variance (e.g.,
TSFEL_0_standard_deviation). The second cluster, in the top right, contains
features that capture different types of autocorrelation structure in the time
series, including linear autocorrelation coefficients (e.g.,
TSFRESH_values_autocorrelation_lag_7 and KATS_seas_acf1). The smaller third
sub-cluster (a sub-group of the autocorrelation cluster) is comprised of two
copies of the exact same feature as described earlier—
SC_FluctAnal_2_rsrangefit_50_1_logi_prop_r1 (which is contained in both
catch22 and tsfeatures). This analysis reveals that the two broad types of
time-series properties that best distinguish the classes of the Bonn EEG
dataset are in fact very simple: variance and linear autocorrelation
structure. This understanding was obtained by comparing the results across six
different open-source feature sets, aided by the ability to inspect the table
of top performing features alongside the clustered feature–feature correlation
plot. However, while the differences between classes in this case were simple,
other, more complex features may perform the strongest on other problems, or
even different pairs of classes within the five-class dataset investigated
here. Identifying when simple features perform well is important as it can
provide interpretable benchmarks for assessing relative performance gains
achieved by more complex and/or less interpretable alternative classifiers.

Figure 5:  A group of variance-sensitive features and a group of
autocorrelation-sensitive features perform the best at distinguishing between
the five classes in the Bonn EEG dataset using the absolute Spearman
correlation coefficient, $|\rho|$, to capture feature–feature similarity. To
aid the identification of similarly performing features, the matrix of
correlation coefficients between features were then organized using
hierarchical clustering (on Euclidean distances with average linkage) along
rows and columns to order the heatmap graphic. This plot was generated by the
code in Listing 8.

Having identified the discriminative features, it is important understand how
they differ amongst the labeled classes of a dataset. This can be achieved by
visualizing the distribution of values for each class for each of the
features. In theft, compute_top_features produces the object ViolinPlots,
where each time series is represented as a point colored by its class label.
Example code is shown in Listing 8, which produces violin plots for all 40 top
features. Here, for visual clarity, we show the violin plots for a selected
feature from the variance-sensitive cluster of features from Fig. 5:
0_standard_deviation from TSFEL (measures the standard deviation); and a
selected feature from the autocorrelation-sensitive cluster of features:
values_autocorrelation_lag_5 from tsfresh (calculates the autocorrelation
coefficient at a time lag of 5 samples). The outputs are shown in Fig. 6.
Consistent with their high classification scores, both features are
informative of class differences. The plot shows that with regards to
autocorrelation structure, we see that ‘eyesClosed’ exhibits the lowest
coefficient at lag 5, while ‘hippocampus’ and ‘epileptogenic’ exhibit the
highest. The plot also highlights that ‘seizure’ time series have increased
standard deviation, consistent with prior work (Fulcher _et al._ , 2013).

Figure 6:  Violin plots (on original feature value scale) of a sample of two
of the top 40 features of all six feature sets in theft for classifying Bonn
EEG groups from the compute_top_features function. Classes differ in their
variance and autocorrelation properties. The features selected for this plot
were determined through the code in Listing 8.

###  2.8 Additional functionality

In addition to the functionality demonstrated here, theft includes a
collection of other functions not demonstrated here, including visualizations
of pairwise correlation matrices (of both feature vectors and raw time-series
values), processing of hctsa-formatted Matlab files, and a number of functions
for investigating and cleaning feature data. These have been omitted from this
article for space, but readers are encouraged to explore them in the detailed
vignette that is included in theft, and in the source code (Henderson and
Bryant, 2022).

###  2.9 Accompanying interactive web application

To provide the functionality of theft to analysts who may not be fluent with R
and who also seek a fast impression of feature values and class
discrimination, we have also developed an accompanying interactive web
application (Henderson, 2022) written in Shiny (Chang _et al._ ,
2020)222https://dynamicsandneuralsystems.shinyapps.io/timeseriesfeaturevis/.
The application allows users to upload a time-series dataset via a drag-and-
drop interface, and the core functionality of theft can then be performed in
the web browser. Most of the graphical and computational functionality
included in theft is presented in the web browser and users can download a
file of the computed time-series features in a tidy format that can be read
into any analysis program. All the graphics available in theft are presented
as interactive graphics by default to further enable an intuitive exploration
of the uploaded dataset. A screenshot of the informative feature
identification page in the web application is displayed in Fig. 7.

Figure 7:  Example screenshot of the interactive web application
implementation of theft. An example of the interactivity of the top feature
identification page is shown. The adjustable parameters on the left are user-
interface renders of the function arguments in the theft package.

##  3 Discussion

Feature-based time-series analysis is a powerful computational tool for
solving problems using sequential (e.g., time-ordered) data. We have
introduced theft, an open-source package for R which implements the
extraction, processing, visualization, and statistical analysis of time-series
features. The value of time-series features stems from their interpretability
and strong connection to theory that can be used to understand empirical
dynamics. theft provides a unified interface to extracting features from six
open-source packages—catch22, feasts, tsfeatures, Kats, tsfresh, and
TSFEL—along with a comprehensive range of analyses to leverage the combined
contributions from all of these packages. For the first time in the free and
open-source software setting, theft provides a full workflow for conducting
feature-based time-series analysis, taking the analyst from feature extraction
through to generating interpretable insights about their data. We demonstrated
theft on the five-class Bonn EEG time-series classification problem (Andrzejak
_et al._ , 2001), in which the full feature-based classification analysis
pipeline—from feature extraction to normalization, classification, and
interpretation of individual features—was achieved using a small number of key
functions in theft. theft can compare feature-set performance and leverage the
combined set of features from all six packages, with in-built techniques like
low-dimensional projections (plot_low_dimension) and feature–feature
correlation matrices (compute_top_features) assisting in interpreting the
patterns detected. Analysts no longer need to construct complex workflows with
multiple software libraries that were not designed to work together—theft
provides a full suite of functionality, but also provides a blueprint for
advanced users to alter and adapt as their research requires.

As new and more powerful features (and feature sets) are developed in the
future, they can be incorporated into theft to enable ongoing assessments of
the types of problems they are best placed to solve. In addition to the
analysis templates provided through functions in theft, there is much
flexibility for users to adapt them or build new functionality for their own
use-cases, such as applying different types of statistical learning algorithms
on extracted feature matrices (e.g., feature selection), or to adapt the
results to different applications such as extrinsic regression (Tan _et al._ ,
2021) or forecasting (Montero-Manso _et al._ , 2020). Future work could also
aim to reduce redundancy from across the combined features towards a new
reduced feature set that combines the most generically informative and unique
features from across the available feature-extraction packages (following the
aims of the catch22 feature set, selected from a library of $>7700$ candidate
features in hctsa Lubba _et al._ (2019)).

##  4 Code Availability

The source code for theft is available on GitHub at
https://github.com/hendersontrent/theft (Henderson and Bryant, 2022).

## References

  * Andrzejak _et al._ (2001) Andrzejak RG, Lehnertz K, Mormann F, Rieke C, David P, Elger CE (2001).  “Indications of Nonlinear Deterministic and Finite-Dimensional Structures in Time Series of Brain Electrical Activity: Dependence on Recording Region and Brain State.”  _Physical Review. E, Statistical, Nonlinear, and Soft Matter Physics_ , 64(6 Pt 1), 061907.  ISSN 1539-3755.  10.1103/PhysRevE.64.061907.
  * Barandas _et al._ (2020) Barandas M, Folgado D, Fernandes L, Santos S, Abreu M, Bota P, Liu H, Schultz T, Gamboa H (2020).  “TSFEL: Time Series Feature Extraction Library.”  _SoftwareX_ , 11, 100456.  ISSN 2352-7110.  10.1016/j.softx.2020.100456.
  * Barbara _et al._ (2022) Barbara NH, Bedding TR, Fulcher BD, Murphy SJ, Van Reeth T (2022).  “Classifying Kepler Light Curves for 12,000 A and F Stars Using Supervised Feature-Based Machine Learning.”  _Monthly Notices of the Royal Astronomical Society_ , p. stac1515.  ISSN 0035-8711.  10.1093/mnras/stac1515.
  * Chang _et al._ (2020) Chang W, Cheng J, Allaire J, Xie Y, McPherson J (2020).  _shiny: Web Application Framework for R_.  R package version 1.5.0, URL https://CRAN.R-project.org/package=shiny.
  * Christ _et al._ (2018) Christ M, Braun N, Neuffer J, Kempa-Liehr AW (2018).  “Time Series FeatuRe Extraction on Basis of Scalable Hypothesis Tests (Tsfresh – A Python Package).”  _Neurocomputing_ , 307, 72–77.  ISSN 0925-2312.  10.1016/j.neucom.2018.03.067.
  * Christ _et al._ (2017) Christ M, Kempa-Liehr AW, Feindt M (2017).  “Distributed and Parallel Time Series Feature Extraction for Industrial Big Data Applications.”  10.48550/arXiv.1610.07717.  1610.07717.
  * Cleveland _et al._ (1990) Cleveland RB, Cleveland WS, McRae JE, Terpenning I (1990).  “STL: A Seasonal-Trend Decomposition Procedure Based on Loess (with Discussion).”  _Journal of Official Statistics_ , 6, 3–73.
  * Day and Edelsbrunner (1984) Day WHE, Edelsbrunner H (1984).  “Efficient Algorithms for Agglomerative Hierarchical Clustering Methods.”  _Journal of Classification_ , 1(1), 7–24.  ISSN 1432-1343.  10.1007/BF01890115.
  * Decat _et al._ (2022) Decat N, Walter J, Koh ZH, Sribanditmongkol P, Fulcher BD, Windt JM, Andrillon T, Tsuchiya N (2022).  “Beyond Traditional Sleep Scoring: Massive Feature Extraction and Data-Driven Clustering of Sleep Time Series.”  _Sleep Medicine_ , 98, 39–52.  ISSN 1389-9457.  10.1016/j.sleep.2022.06.013.
  * Facebook Infrastructure Data Science (2021) Facebook Infrastructure Data Science (2021).  “Kats.”  URL https://facebookresearch.github.io/Kats/.
  * Fulcher (2018) Fulcher BD (2018).  “Feature-Based Time-Series Analysis.”  In _Feature Engineering for Machine Learning and Data Analytics_. CRC Press.  ISBN 978-1-315-18108-0.
  * Fulcher _et al._ (2012) Fulcher BD, Georgieva AE, Redman CWG, Jones NS (2012).  “Highly Comparative Fetal Heart Rate Analysis.”  In _2012 Annual International Conference of the IEEE Engineering in Medicine and Biology Society_ , pp. 3135–3138.  ISSN 1558-4615.  10.1109/EMBC.2012.6346629.
  * Fulcher and Jones (2014) Fulcher BD, Jones NS (2014).  “Highly Comparative Feature-Based Time-Series Classification.”  _IEEE Transactions on Knowledge and Data Engineering_ , 26(12), 3026–3037.  ISSN 1041-4347, 1558-2191, 2326-3865.  10.1109/TKDE.2014.2316504.  1401.3531.
  * Fulcher and Jones (2017) Fulcher BD, Jones NS (2017).  “Hctsa: A Computational Framework for Automated Time-Series Phenotyping Using Massive Feature Extraction.”  _Cell Systems_ , 5(5), 527–531.e3.  ISSN 2405-4712.  10.1016/j.cels.2017.10.001.
  * Fulcher _et al._ (2013) Fulcher BD, Little MA, Jones NS (2013).  “Highly Comparative Time-Series Analysis: The Empirical Structure of Time Series and Their Methods.”  _Journal of The Royal Society Interface_ , 10(83), 20130048\.  10.1098/rsif.2013.0048.
  * Fulcher _et al._ (2020) Fulcher BD, Lubba CH, Sethi SS, Jones NS (2020).  “A Self-Organizing, Living Library of Time-Series Data.”  _Scientific Data_ , 7(1), 213.  ISSN 2052-4463.  10.1038/s41597-020-0553-0.
  * Harris (2021) Harris BJ (2021).  _Catch22.jl_.  https://doi.org/10.5281/zenodo.5030712.  V0.2.1.
  * Henderson (2021) Henderson T (2021).  _Rcatch22: Calculation of 22 CAnonical Time-Series CHaracteristics_.  R package version 0.1.12.
  * Henderson (2022) Henderson T (2022).  “hendersontrent/theft-webtool: v0.1.1.”  10.5281/ZENODO.6656286.  URL https://zenodo.org/record/6656286.
  * Henderson and Bryant (2022) Henderson T, Bryant AG (2022).  “hendersontrent/theft: v0.3.9.7.”  10.5281/ZENODO.6650876.  URL https://zenodo.org/record/6650876.
  * Henderson and Fulcher (2021) Henderson T, Fulcher BD (2021).  “An Empirical Evaluation of Time-Series Feature Sets.”  In _2021 International Conference on Data Mining Workshops (ICDMW)_ , pp. 1032–1038.  ISSN 2375-9259.  10.1109/ICDMW53433.2021.00134.
  * Hyndman _et al._ (2020) Hyndman R, Kang Y, Montero-Manso P, Talagala T, Wang E, Yang Y, O’Hara-Wild M (2020).  _tsfeatures: Time Series Feature Extraction_.  R package version 1.0.2, URL https://CRAN.R-project.org/package=tsfeatures.
  * Jolliffe (2002) Jolliffe IT (2002).  _Principal Component Analysis_.  Springer Series in Statistics. Springer-Verlag, New York.  ISBN 978-0-387-95442-4.  10.1007/b98835.
  * Kao _et al._ (2021) Kao LJ, Chiu CC, Wang HJ, Ko CY (2021).  “Prediction of Remaining Time on Site for E-Commerce Users: A SOM and Long Short-Term Memory Study.”  _Journal of Forecasting_ , n/a(n/a).  ISSN 1099-131X.  10.1002/for.2771.
  * Kodra _et al._ (2011) Kodra E, Chatterjee S, Ganguly AR (2011).  “Exploring Granger Causality between Global Average Observed Time Series of Carbon Dioxide and Temperature.”  _Theoretical and Applied Climatology_ , 104(3), 325–335.  ISSN 1434-4483.  10.1007/s00704-010-0342-3.
  * Kuhn (2020) Kuhn M (2020).  _caret: Classification and Regression Training_.  R package version 6.0-86, URL https://CRAN.R-project.org/package=caret.
  * Liu _et al._ (2020) Liu G, Li L, Zhang L, Li Q, Law SS (2020).  “Sensor Faults Classification for SHM Systems Using Deep Learning-Based Method with Tsfresh Features.”  _Smart Materials and Structures_ , 29(7), 075005.  ISSN 0964-1726.  10.1088/1361-665X/ab85a6.
  * Lubba _et al._ (2019) Lubba CH, Sethi SS, Knaute P, Schultz SR, Fulcher BD, Jones NS (2019).  “Catch22: CAnonical Time-series CHaracteristics.”  _Data Mining and Knowledge Discovery_ , 33(6), 1821–1852.  ISSN 1573-756X.  10.1007/s10618-019-00647-x.
  * Markicevic _et al._ (2020) Markicevic M, Fulcher BD, Lewis C, Helmchen F, Rudin M, Zerbi V, Wenderoth N (2020).  “Cortical Excitation:Inhibition Imbalance Causes Abnormal Brain Network Dynamics as Observed in Neurodevelopmental Disorders.”  _Cerebral Cortex_ , 30(9), 4922–4937.  ISSN 1047-3211.  10.1093/cercor/bhaa084.
  * Montero-Manso _et al._ (2020) Montero-Manso P, Athanasopoulos G, Hyndman RJ, Talagala TS (2020).  “FFORMA: Feature-based Forecast Model Averaging.”  _International Journal of Forecasting_ , 36(1), 86–92.  ISSN 0169-2070.  10.1016/j.ijforecast.2019.02.011.
  * Ojala and Garriga (2009) Ojala M, Garriga GC (2009).  “Permutation Tests for Studying Classifier Performance.”  In _2009 Ninth IEEE International Conference on Data Mining_ , pp. 908–913. IEEE, Miami Beach, FL, USA.  ISBN 978-1-4244-5242-2.  10.1109/ICDM.2009.108.
  * Paul _et al._ (2021) Paul A, McLendon H, Rally V, Sakata JT, Woolley SC (2021).  “Behavioral Discrimination and Time-Series Phenotyping of Birdsong Performance.”  _PLOS Computational Biology_ , 17(4), e1008820.  ISSN 1553-7358.  10.1371/journal.pcbi.1008820.
  * Sorzano _et al._ (2014) Sorzano COS, Vargas J, Montano AP (2014).  “A Survey of Dimensionality Reduction Techniques.”  10.48550/arXiv.1403.2877.  1403.2877.
  * Subasi and Ismail Gursoy (2010) Subasi A, Ismail Gursoy M (2010).  “EEG Signal Classification Using PCA, ICA, LDA and Support Vector Machines.”  _Expert Systems with Applications_ , 37(12), 8659–8666.  ISSN 0957-4174.  10.1016/j.eswa.2010.06.065.
  * Tan _et al._ (2021) Tan CW, Bergmeir C, Petitjean F, Webb GI (2021).  “Time Series Extrinsic Regression.”  _Data Mining and Knowledge Discovery_ , 35(3), 1032–1060.  ISSN 1573-756X.  10.1007/s10618-021-00745-9.
  * Van Der Donckt _et al._ (2022) Van Der Donckt J, Van Der Donckt J, Deprost E, Van Hoecke S (2022).  “Tsflex: Flexible Time Series Processing & Feature Extraction.”  _SoftwareX_ , 17, 100971.  ISSN 2352-7110.  10.1016/j.softx.2021.100971.
  * van der Maaten and Hinton (2008) van der Maaten L, Hinton G (2008).  “Visualizing Data Using T-SNE.”  _Journal of Machine Learning Research_ , 9(86), 2579–2605.  ISSN 1533-7928.
  * West _et al._ (1999) West M, Prado R, Krystal AD (1999).  “Evaluation and Comparison of EEG Traces: Latent Structure in Nonstationary Time Series.”  _Journal of the American Statistical Association_ , 94(446), 375–387.  ISSN 0162-1459.  10.1080/01621459.1999.10474128.
  * Wickham (2014) Wickham H (2014).  “Tidy Data.”  _Journal of Statistical Software_ , 59(1), 1–23.  ISSN 1548-7660.  10.18637/jss.v059.i10.
  * Wickham _et al._ (2019) Wickham H, Averick M, Bryan J, Chang W, McGowan LD, François R, Grolemund G, Hayes A, Henry L, Hester J, Kuhn M, Pedersen TL, Miller E, Bache SM, Müller K, Ooms J, Robinson D, Seidel DP, Spinu V, Takahashi K, Vaughan D, Wilke C, Woo K, Yutani H (2019).  “Welcome to the Tidyverse.”  _Journal of Open Source Software_ , 4(43), 1686.  ISSN 2475-9066.  10.21105/joss.01686.
  * Yang _et al._ (2021) Yang Z, Abbasi IA, Mustafa EE, Ali S, Zhang M (2021).  “An Anomaly Detection Algorithm Selection Service for IoT Stream Data Based on Tsfresh Tool and Genetic Algorithm.”  _Security and Communication Networks_ , 2021, 6677027.  ISSN 1939-0114.  10.1155/2021/6677027.