\section{Introduction}
\label{sec:introduction}

\IEEEPARstart{I}{nterval-valued} data have gained increasing attention across a variety of contexts, as they offer a systematic way to capture information complete with an intrinsic representation of range or uncertainty in each individual `measurement', which is not possible using point-values, such as numbers or ranks. Such data may arise due to imprecision and uncertainty in measurement (e.g., sensor data), inherent uncertainty of outcome (e.g., estimations of stock prices ~\cite{lin2020symbolic}, climate/weather forecasting~\cite{sun2015normal}), or inherent vagueness or nuance (e.g., in linguistic terms~\cite{liu2008encoding}). A growing body of evidence also suggests that subjective judgement of humans (e.g., experts, consumers) may be better represented by intervals~\cite{wagner2015interval,Ellerby2020CapturingIU}, where the interval width/size captures the response range as a conjunctive set (e.g., the reals between 2 and 4), or degree of uncertainty in respect to an estimate or rating ---a disjunctive set (e.g., confidence interval) \cite{kabir2020similarity}.

Regression analysis for estimating one or more variables is an important task in data analysis, pattern recognition, and machine learning under the wider umbrella of AI. While regression is most commonly applied to numeric data, researchers have recently begun to develop methods of applying (linear) regression analysis to interval-valued data~\cite{brito2007modelling}. To date, existing models have been developed either within the framework of random sets or symbolic data analysis (SDA)~\cite{yeh2017tree,gonzalez2013constrained}. In the random set framework, intervals are viewed as compact convex sets in $\mathbb{R}$, and their interaction is modelled according to set arithmetic~\cite{gonzalez2007least,blanco2011estimation,boukezzoula2011midpoint,blanco2013set,garcia2020multiple}. On the other hand, in the SDA setting, the classical regression model is expanded to interval-valued data where intervals are treated as bivariate vectors~\cite{billard2000regression,billard2002symbolic,neto2008centre,neto2010constrained,wang2012linear,sun2015linear,souza2017parametrized,hao2017constrained}.

Current linear regression models for intervals use different \emph{reference points}, such as the regressors' center values, their lower and upper bounds, or both their center and range (width), to estimate the interval-valued regressand~\cite{souza2017parametrized}. All earlier interval regression approaches~\cite{billard2000regression,billard2002symbolic,neto2008centre} risk `flipping' the interval bounds or generating a negative interval range, and thus breaking mathematical coherence (i.e., the value of the left endpoint is greater than the value of the right endpoint, violating the definition of the interval). The more recent approaches~\cite{neto2010constrained,wang2012linear,sun2015linear,souza2017parametrized} impose either positivity restrictions on parameters or design the regression model so as to ensure that the lower bound does not exceed the upper bound. For instance, Neto and Carvalho~\cite{neto2010constrained} apply the Lawson and Hanson's algorithm (LHA)~\cite{lawson1974},  Wang et al.~\cite{wang2012linear} use Moore’s linear combination~\cite{moore1966interval}, and  Souza et al.~\cite{souza2017parametrized} the Box-Cox transformation~\cite{box1964analysis} to guarantee mathematical coherence. Although Sun and Ralescu~\cite{sun2015linear} consider positivity restrictions on parameters in their model setting, they do not provide a method to maintain these restrictions. A detailed review and analysis of these models is provided in Section~\ref{sec:background} and the issue is discussed throughout the paper.

An important aspect is that most current approaches, including the above, restrict parameters based on the assumption that an increase in uncertainty\footnote{Throughout the remainder of this paper, for simplicity, we refer to the range of the interval as representing uncertainty.} of a regressor (i.e., input) leads to greater uncertainty in the dependent variable or regressand; however, there are plausible cases in the real world where this positive correlation does not hold. For example, in cyber-security, uncertainty around the vulnerability of a system may lead to certainty in the need to invest resources into controls and thus mitigating the risk of the system being compromised, that is, a negative correlation. Similarly, in a multivariate case, consider an economics example, where greater uncertainty in future product prices and interest rates may cause more certainty in respect to present consumption choices.

For illustration, Fig.~\ref{fig1} presents a synthetic example of two interval-valued data sets where in Set-1, the range (uncertainty) of the regressand, $Y$ increases with the range of regressor, $X$---the more common case as outlined above, e.g., price of cars vs their horsepower. In contrast, for Set-2, the range of regressand, $Y$ increases irrespective of the range of regressor, $X$---a less common case, e.g., the wealth of people in respect to their age. Both examples will be discussed in detail in Section~\ref{subsec:syn_results}.

\begin{figure}
\begin{subfigure}{0.24\textwidth}
\includegraphics[width=0.95\textwidth]{figure1a.png}
\caption{Set-1}
\end{subfigure}
\begin{subfigure}{0.24\textwidth}
\includegraphics[width=0.95\textwidth]{figure1b.png}
\caption{Set-2}
\end{subfigure}
\caption{2D visualization of two different interval-valued sets.}
\label{fig1}
\end{figure}

Aspects such as the above highlight the complexity but also the potential of interval-valued data, in directly measuring uncertainty-at-source and making it accessible for reasoning. With intervals encoding additional information by comparison with numbers, the choice of a `universally best' analysis approach, as can often be found for numeric cases, is highly challenging. Instead, the selection of appropriate techniques for a given interval-valued dataset---based on its properties---is essential. Ferson~et~al., in~\cite{ferson2007experimental}, discuss the complexity of analysing interval-valued data in detail, elucidating key properties for datasets, such as whether the intervals exhibit strong overlap (i.e., `puffy'~\cite{ferson2007experimental}: wider, more uncertain intervals) or very little overlap (i.e., `skinny'~\cite{ferson2007experimental}: narrow, less uncertain); whether there are outliers; whether intervals are highly scattered, densely placed or even nested within each other. Beyond these properties, more fundamental questions arise from the semantic meaning of the given intervals, as already alluded to above; whether they represent an epistemic, disjunctive set, with the interval describing the range within which a correct answer is contained (cf. confidence intervals), or indeed an ontic, conjunctive set with a true range, such as the real numbers between 1 and 5, or a~time span~\cite{couso2014statistical}. For conciseness, we do not discuss the latter further in this manuscript, but will revisit it, and its impact on the analysis of AI and machine learning more generally, in future publications.

To support the broader adoption of intervals as a fundamental data-type in AI, this paper puts forward a detailed analysis of existing linear regression methods in respect to interval-valued datasets with different properties. Succinct illustrations of the relevant properties are provided through synthetic examples, including based on canonical examples inspired by the `puffy' and `skinny' data sets from~\cite{ferson2007experimental}, before addressing a set of real-world interval-valued data. The analyses focus on how closely different linear regression methods fit interval-valued datasets with the given sets of properties. Further, we discuss how the approaches do or do not ensure mathematical coherence, as alluded to above, and the impacts of the adopted approaches on ensuring coherence of model-fit-to-data. In addition, we propose a weaker set of restrictions for the Linear Model (LM)~\cite{sun2015linear}, enhancing its general suitability for real-world data without compromising mathematical coherence.
Specifically, we propose a selective approach to applying restrictions, letting the model's parameters vary freely unless the mathematical coherence is being compromised.

Based on the evaluation of all methods, in terms of both model-fit and mathematical coherence, using both synthetic and real-world data sets, the paper puts forward practical recommendations for linear regression applied to interval-valued data. It also way-points future research on leveraging interval-valued data in AI more generally.

Succinctly, the contributions include the following:
\begin{itemize}
\item Introduction of interval-valued data as a powerful data-type in AI, including detailed discussion of both potential benefits and challenges.
\item Full review and comparative analysis of state-of-the-art approaches to interval-valued regression, for synthetic and real-world data sets.
\item Advancement of the notion of realistic regression behaviours, highlighting the real-world complexity of actual relationships between interval regressor(s) and regressand(s).
\item Extension of the Linear Model~\cite{sun2015linear}, by allowing its parameters to vary freely unless mathematical coherence breaks down, and thus enhancing its suitability for real-world applications.
\item Introduction of novel intuitive visualisation approaches making interval regression results tangible and accessible, vital for their use by decision makers and data scientists.
\item Recommendations on which regression techniques are the most suitable and when.
\item Charting of the gaps in current regression and more general AI techniques that inhibit wider adoption of interval-valued data, and which offer a rich area for future high-impact research.
\end{itemize}

The paper is organized as follows: Section~\ref{sec:background} reviews linear regression methods based on vector/matrix representations of intervals in the regression process. Section~\ref{sec:methodology} discusses adaptations of the LM method~\cite{sun2015linear} for real-world applications, followed by interval-valued data sets (synthetic and real-world) and evaluation metrics. Section~\ref{sec:results} demonstrates the estimation performance of these regression models in respect to datasets. Section~\ref{sec:discussion} discusses suitability of the existing regression approaches to various features of interval data. Lastly,  Section~\ref{sec:conclusion} concludes the paper with future works. Table~\ref{tab:acronym} presents a list of acronyms and notation used in this paper to assist the reader.

\begin{table}
\centering
\caption{Acronyms and Notation}
\begin{tabular}{ll} \hline\hline
CM & Center Method \cite{billard2000regression}\\
MinMax & MinMax Method \cite{billard2002symbolic}\\
CRM & Center and Range Method \cite{neto2008centre}\\
CCRM & Constrained Center and Range Method \cite{neto2010constrained}\\
CIM & Complete Information Method \cite{wang2012linear}\\
LM & Linear Model \cite{sun2015linear}\\
PM & Parametrized Model \cite{souza2017parametrized}\\
LM$_c$ & Constrained Linear Model \\
LM$_w$ & Weakly Constrained Linear Model\vspace{0.5mm}\\\hline
$\overline{a}$ & Interval $\{\overline{a}\subseteq\mathbb{R}:\overline{a}=[a^-,a^+],~a^-\leq a^+\}$ \\
$a^w$ & Range of Interval, $a^w=|\overline{a}|=|a^+-a^-|$ \\
$a^c$ & Center of Interval, $a^c=\frac{(a^++a^-)}{2}$  \\
$a^r$ & Radius of Interval $a^r=\frac{a^w}{2}$ (Half Range)\\
$Y$ & Interval-valued Regressand\\
$\hat{Y}$ & Estimated $Y$\\
$X$ & Interval-valued Regressor\\
$\mu_{Y^w}$ and $\sigma_{Y^w}$ & Mean and Standard Deviation of Range of $Y$\\
$\mu_{X^w}$ and $\sigma_{X^w}$ & Mean and Standard Deviation of Range of $X$\\
$\mu_{Y^c}$ and $\sigma_{Y^c}$ & Mean and Standard Deviation of Center of $Y$\\
$\mu_{X^c}$ and $\sigma_{X^c}$& Mean and Standard Deviation of Center of $X$\\
\hline
\end{tabular}
\label{tab:acronym}
\end{table}

\section{Background}
\label{sec:background}
 A closed interval $\overline{a}$ is characterized by its two endpoints, $a^-$ and $a^+$ with $a^-\leq a^+$ where $a^-$ and $a^+$ are respectively the
lower and upper bounds of $\overline{a}$~\cite{Beliakov2014penalty}. Generally, $\overline{a}$ is denoted as $\overline{a}=[a^-,a^+]$, and its range is given by  $\left|\overline{a}\right|=a^w=\left|a^+-a^-\right|$. Alternatively, the same interval can be represented as $\overline{a}=[a^c-a^r,a^c+a^r]$ with its center, $a^c=\frac{(a^+ + a^-)}{2}$ and radius, $a^r=\frac{a^w}{2}$. In this section, we review well-known linear regression methods for interval-valued data where they are treated as bivariate vectors in the regression process. Here, $Y=\{\overline{y}_1,\overline{y}_2,\dots,\overline{y}_n\}$ is an interval-valued regressand (dependent variable) with $n$ observations where $\overline{y}_i=[y_i^-,y_i^+]$, $1\leq i\leq n$. Further, $\{X_1,X_2,\dots,X_p\}$ presents $p$ ($\geq1$) interval-valued regressors (independent variables) where each $X_j$ also has $n$ observations, $X_j=\{\overline{x}_{j1},\overline{x}_{j2},\dots,\overline{x}_{jn}\}$ with $\overline{x}_{ji} = [x_{ji}^-,x_{ji}^+]$, $1\leq i\leq n$, $1\leq j\leq p$. In addition, the estimated value of $Y$ is given by  $\hat{Y}=\{\hat{\overline{y}}_1,\hat{\overline{y}}_2,\dots,\hat{\overline{y}}_n\}$.

\subsection{The Center Method (2000)}
\label{CM_subsection}
Billard and Diday~\cite{billard2000regression} proposed a linear regression model, known as the Center Method (CM) in 2000, to fit the center of both regressand and regressors, and then apply this model to the lower and upper bounds of regressors to separately estimate the upper and lower bounds of the regressand. Equation~\eqref{CM_lower_upper} defines the CM regression model,
\begin{equation}\label{CM_lower_upper}
\begin{split}
    y_i^-=\beta_0^c+\sum_{j=1}^p\beta_j^c x_{ji}^-+\epsilon_i^-\\
    y_i^+=\beta_0^c+\sum_{j=1}^p\beta_j^c x_{ji}^++\epsilon_i^+
\end{split}
\end{equation}
where $y_i^-$ and $y_i^+$ are the lower and upper bounds of $\overline{y}_i$. Similarly,  $x_{ji}^-$ and $ x_{ji}^+$ are the the lower and upper bounds of $\overline{x}_{ji}$. $\beta_j^c\in \mathbb{R}$ are the regression coefficients for the center estimation, $0\leq j\leq p$ while $\epsilon_i^-$ and $\epsilon_i^+$ are error terms. Both lower and upper bound models of~\eqref{CM_lower_upper} use the same coefficients, therefore, they can be captured together using~\eqref{CM_center}.
\begin{equation}\label{CM_center}
    y_i^c=\beta_0^c+\sum_{j=1}^p\beta_j^c x_{ji}^c+\epsilon_i^c,
\end{equation}
where $y_i^c$ and $x_{ji}^c$ are the center of $\overline{y}_i$ and $\overline{x}_{ji}$ respectively and $\epsilon_i^c=(\epsilon_i^-+\epsilon_i^+)/2$ is error term. Equation~\eqref{CM_center} can also be expressed for the total set of $n$ observations in the matrix format as $Y^c=X^{c}\beta^c+\epsilon^c$, where

\resizebox{.94\columnwidth}{!}{
\begin{minipage}{0.9\columnwidth}
\begin{align*}\label{CM_matrix_form}
    \begin{split}
     Y^c=(y_1^c\quad y_2^c\quad \dots\quad y_n^c)^T,\quad \beta^c=(\beta_0^c\quad \beta_1^c\quad\dots\quad\beta_p^c)^T\\
     \epsilon^c=(\epsilon_1^c\quad \epsilon_2^c\quad\dots\quad\epsilon_n^c)^T, \text{ and } X^c=\begin{pmatrix}
                     1 & x_{11}^c & \dots & x_{p1}^c\\
                     1 & x_{12}^c & \dots & x_{p2}^c\\
                     \vdots & \vdots & \dots &\vdots\\
                     1 & x_{1n}^c &\dots & x_{pn}^c\\
               \end{pmatrix}.
\end{split}
\end{align*}
\end{minipage}}

The coefficients, $\beta^c$ are estimated using the least squares (LS) method~\cite{draper1998applied}, defined in~\eqref{eq:Beta_CM}.
\begin{equation}\label{eq:Beta_CM}
    \hat{\beta}^c=((X^c)^T X^c)^{-1}(X^c)^TY^c
\end{equation}
By using $\hat{\beta}^c$, ~\eqref{Estimated_CM_lower_upper} estimates the bounds of regressand, $Y$.
\begin{equation}\label{Estimated_CM_lower_upper}
\begin{split}
    \hat{y}_i^-=\hat{\beta}_0^c+\sum_{j=1}^p\hat{\beta}_j^c x_{ji}^-\\
    \hat{y}_i^+=\hat{\beta}_0^c+\sum_{j=1}^p\hat{\beta}_j^c x_{ji}^+.
\end{split}
\end{equation}

The CM model is simple, as it follows a standard regression approach applied to the center of the intervals and then the resulting model is applied both to upper and  lower endpoints. However, by doing so, it risks and often fails to maintain mathematical coherence \cite{neto2008centre}, that is, the estimated lower bound of $\overline{y}_i$ becomes larger than the estimated upper bound.

\subsection{The MinMax Method (2002)} \label{MM_subsection}
Billard and Diday~\cite{billard2002symbolic} presented the MinMax method in 2002 which directly uses the lower and upper bounds of the regressors to separately estimate the lower and upper bounds of the regressand. Equation~\eqref{MinMax} defines the MinMax regression model for the lower and upper bounds of the regressand,
\begin{equation}\label{MinMax}
\begin{split}
    y_i^-=\beta_0^-+\sum_{j=1}^p\beta_j^- x_{ji}^-+\epsilon_i^-\\
    y_i^+=\beta_0^++\sum_{j=1}^p\beta_j^+ x_{ji}^++\epsilon_i^+
\end{split}
\end{equation}
where $\beta_j^-$ and $\beta_j^+$ are regression coefficients for the lower and upper bounds respectively ($\beta_j^-,\beta_j^+\in \mathbb{R}$, $0\leq j\leq p$). In matrix notation, the lower bound model of~\eqref{MinMax} with all $n$ observations is expressed as $Y^-=X^{-}\beta^-+\epsilon^-$, where

\resizebox{.94\columnwidth}{!}{
  \begin{minipage}{0.9\columnwidth}
 \begin{align*}\label{MinMax_lower_matrix_form}
    \begin{split}
      Y^-=(y_1^-\quad y_2^-\quad \dots \quad y_n^-)^T,\quad \beta^-=(\beta_0^-\quad \beta_1^-\quad \dots \quad\beta_p^-)^T,\\
      \epsilon^-=(\epsilon_1^-\quad \epsilon_2^-\quad \dots \quad\epsilon_n^-)^T, \text{ and } X^-=\begin{pmatrix}
                     1 & x_{11}^- & \dots & x_{p1}^-\\
                      1 & x_{12}^- & \dots & x_{p2}^-\\
                     \vdots & \vdots & \dots &\vdots\\
                      1 & x_{1n}^- &\dots & x_{pn}^-\\
               \end{pmatrix}.
 \end{split}
 \end{align*}
 \end{minipage}
 }

Like the CM method, \eqref{eq:Beta_MinMax_lower} estimates the coefficient ($\beta^-$) for the lower bound model of~\eqref{MinMax} using the LS method~\cite{draper1998applied}.
\begin{equation}\label{eq:Beta_MinMax_lower}
    \hat{\beta}^-=((X^-)^T X^-)^{-1}(X^-)^TY^-
\end{equation}
Now, using the estimated $\hat{\beta}^-$, \eqref{Estimated_MinMax_lower} estimates the lower bounds for the regressand, $Y$.
\begin{equation}\label{Estimated_MinMax_lower}
    \hat{y}_i^-=\hat{\beta}_0^-+\sum_{j=1}^p\hat{\beta}_j^- x_{ji}^-
\end{equation}

Similarly, the upper bound model of \eqref{MinMax} is expressed as $Y^+=X^{+}\beta^{+}+\epsilon^+$, where

 \resizebox{0.94\columnwidth}{!}{
   \begin{minipage}{0.9\columnwidth}
 \begin{align*
     \begin{split}
     Y^+=(y_1^+\quad y_2^+\quad\dots\quad y_n^+)^T,\quad \beta^+=(\beta_0^+\quad \beta_1^+\quad\dots\quad\beta_p^+)^T,\\
      \epsilon^+=(\epsilon_1^+\quad \epsilon_2^+\quad\dots\quad\epsilon_n^+)^T, \text{ and } X^+=\begin{pmatrix}
                      1 & x_{11}^+ & \dots & x_{p1}^+\\
                      1 & x_{12}^+ & \dots & x_{p2}^+\\
                      \vdots & \vdots & \dots &\vdots\\
                      1 & x_{1n}^+ &\dots & x_{pn}^+\\
               \end{pmatrix}.
 \end{split}
 \end{align*}
 \end{minipage}
 }

Equation~\eqref{eq:Beta_MinMax_upper} estimates the coefficients for the upper bound of \eqref{MinMax} by the LS method~\cite{draper1998applied}.
\begin{equation}\label{eq:Beta_MinMax_upper}
    \hat{\beta}^+=((X^+)^T X^+)^{-1}(X^+)^TY^+
\end{equation}
Now, using $\hat{\beta}^+$, \eqref{Estimated_MinMax_upper} estimates the upper bounds of the regressand, $Y$.
\begin{equation}\label{Estimated_MinMax_upper}
      \hat{y}_i^+=\hat{\beta}_0^++\sum_{j=1}^p\hat{\beta}_j^+ x_{ji}^+
\end{equation}

In the MinMax method, the use of two separate models to estimate the lower and upper bounds of the regressand improves the model fitness and interpretation compared to the CM method~\cite{neto2008centre}. However, it still does not guarantee mathematical coherence on the estimated bounds~\cite{neto2008centre}. Further, its estimation performance can be reduced if there is no dependency between the bounds of regressand and regressors~\cite{souza2017parametrized}.

\subsection{The Center and Range Method (2008)}
\label{CRM_subsection}
Neto and Carvalho \cite{neto2008centre} extended the CM method \cite{billard2000regression} in 2008 to the Center and Range Method (CRM), by considering both interval center and half the `range' (i.e., the radius) of the regressor and regressand variables. They build two separate regression models---one for the centers and other for the ranges (or radii). Equation~\eqref{CRM_center} is used to model the interval centers.
\begin{equation}\label{CRM_center}
    y_i^c=\beta_0^c+\sum_{j=1}^p\beta_j^c x_{ji}^c+\epsilon_i^c
\end{equation}
where $y_i^c$ and $x_{ji}^c$ are the centers $\overline{y}_i$ and $\overline{x}_{ji}$ respectively and $\epsilon_i^c$ is error term. The CRM method follows the same matrix format of the CM method to represent the center model of~\eqref{CRM_center} and applies the same strategy of the CM method to estimate the coefficients $\beta_j^c$ (see~\eqref{eq:Beta_CM}). It then estimates the center values for the regressand, $Y$ using the estimated $\hat{\beta}_j^c$ by~\eqref{Estimated_CRM_center}.
\begin{equation}\label{Estimated_CRM_center}
      \hat{y}_i^c=\hat{\beta}_0^c+\sum_{j=1}^p\hat{\beta}_j^c x_{ji}^c
\end{equation}

Similarly, the interval ranges are modelled by the CRM method using~\eqref{CRM_range}, \begin{equation}\label{CRM_range}
    y_i^r=\beta_0^r+\sum_{j=1}^p\beta_j^r x_{ji}^r+\epsilon_i^r,
\end{equation}
where $y_i^r$ and $x_{ji}^r$ are the radius of $\overline{y}_i$ and $\overline{x}_{ji}$ respectively. $\beta_j^r\in\mathbb{R}$ are the coefficients for the range estimations, $0\leq j\leq p$ and $\epsilon_i^r$ is error term. Again, following the same matrix notation, the coefficients, $\beta_j^r$ are estimated by~\eqref{eq:Beta_CRM_range},
\begin{equation}\label{eq:Beta_CRM_range}
    \hat{\beta}^r=((X^r)^T X^r)^{-1}(X^r)^TY^r
\end{equation}
and using the estimated $\hat{\beta}^r$, \eqref{Estimated_CRM_range} computes the ranges of the regressand, $Y$.
\begin{equation}\label{Estimated_CRM_range}
      \hat{y}_i^r=\hat{\beta}_0^r+\sum_{j=1}^p\hat{\beta}_j^r x_{ji}^r
\end{equation}

Finally, the CRM method estimates the bounds of the regressand from the estimated $\hat{y}_i^c$ and $\hat{y}_i^r$ using~\eqref{CRM_lower_upper}.
\begin{equation}\label{CRM_lower_upper}
    \hat{y}_i^-=\hat{y}_i^c-\hat{y}_i^r \text{  and  }  \hat{y}_i^+=\hat{y}_i^c+\hat{y}_i^r.
\end{equation}

The CRM model provides better estimation when there is a linear dependency between the ranges of regressand and regressors~\cite{neto2010constrained,souza2017parametrized}. If this is not the case, it does not ensure mathematical coherence as~\eqref{Estimated_CRM_range} can produce negative $y_i^r$~\cite{souza2017parametrized}.

\subsection{The Constrained Center and Range Method (2010)} \label{CCRM_subsection}
Lima Neto and Carvalho~\cite{neto2010constrained} later refined the CRM model~\cite{neto2008centre} in 2010 to ensure mathematical coherence. This approach is known as the Constrained Center and Range Method (CCRM) where a positivity restriction is enforced on the coefficients which are estimated in respect to the relationship of the radii of regressand and regressor variables. In other words, the overall estimation process in the CCRM method remains the same as the CRM model (see \eqref{CRM_center} to \eqref{CRM_lower_upper} in Section~\ref{CRM_subsection}) with positive constraints on the range coefficients, $\beta_j^r$. The CCRM model uses an iterative algorithm proposed by Lawson and Hanson in~\cite{lawson1974} to ensure $\beta_j^r\geq0$. For more details, please see~\cite{neto2010constrained}. Finally, it estimates the lower and upper bounds of the regressand using~\eqref{CRM_lower_upper}.

The positivity restriction on the range coefficients in the CCRM model guarantees mathematical coherence of the estimated bounds~\cite{neto2010constrained}. However, it can lead to biased estimation outcomes, i.e., a poor fitness to actual intervals~\cite{neto2010constrained}. Hence, it is recommended by~\cite{neto2010constrained} to apply the CCRM model only when the CRM model fails to maintain mathematical coherence.

\subsection{The Complete Information Method (2012)}
\label{CIM_subsection}
Wang et al.~\cite{wang2012linear} presented the Complete Information Method (CIM) in 2012 with a focus on considering all the information contained within the intervals. The CIM model considers each interval observation of regressand and regressor variables as a hyper-cube (if a single-valued regressor and associated regressand are considered, each observation effectively reflects a rectangular, as visualised in Fig.~\ref{fig:synthetic_dataset_A}), and the regression model is built on these hyper-cubes. Further, it adopts Moore’s linear combination algorithm~\cite{moore1966interval} to ensure the mathematical consistency of the bounds. In this model, $Y$ is presented by a linear combination of $X_{j}$,
\begin{equation}\label{CIM}
Y=\beta_0 1_n+\sum_{j=1}^p\beta_j X_{j} +\epsilon_i,
\end{equation}
where $\beta_j$ represents the coefficients with $0\leq j\leq p$ and $1_n$ is an n-dimensional column constant vector of ones. Equation~\eqref{CIM_lower_upper} defines the lower and upper bound estimations of $Y$.
\begin{equation}\label{CIM_lower_upper}
\begin{split}
    y_i^-=\beta_0+\sum_{j=1}^p\beta_j\left(\tau_j x_{ji}^-+(1-\tau_j)x_{ji}^+\right)\\
    y_i^+=\beta_0+\sum_{j=1}^p\beta_j\left((1-\tau_j) x_{ji}^-+\tau_j x_{ji}^+\right)
\end{split}
\end{equation}
where,
\[\tau_j = \left\{
  \begin{array}{lr}
    0, & \text{if } \beta_j\leq 0\\
    1, & \text{otherwise. }
  \end{array}\right.
\]
Using~\eqref{CIM}, the LS method~\cite{draper1998applied} of~\eqref{CIM_lower_upper} generates the following linear system,
\begin{equation}\label{CIM1}
    M\beta = b
\end{equation}
where,
\begin{equation*}
M=\begin{pmatrix}
         \langle 1,1\rangle & \langle 1,X_1\rangle& \dots & \langle 1,X_p\rangle\\
         \langle X_1,1\rangle & \langle X_1,X_1\rangle& \dots & \langle X_1,X_p\rangle\\
          \vdots & \vdots & \dots & \vdots\\
         \langle X_p,1\rangle& \langle X_p,1\rangle&\dots & \langle X_p,X_p\rangle\\
   \end{pmatrix},
\end{equation*}
\begin{equation*}
 \beta=\begin{pmatrix}
         \beta_0\\
         \beta_1\\
         \vdots\\
         \beta_p\\
       \end{pmatrix},
 \text{ and }
 b=\begin{pmatrix}
      \langle 1,Y\rangle\\
      \langle X_1,Y\rangle\\
      \vdots\\
      \langle X_p,Y\rangle\\
   \end{pmatrix}.
\end{equation*}
\noindent Here, $\langle X_a,X_b\rangle$ is the inner product between $X_a$ and $X_b$ ($0\leq a,b\leq p$ and $X_0=1$)  which is computed using~\eqref{CIM_More}.
\begin{equation}\label{CIM_More}
\resizebox{.48\textwidth}{!}{$
 \langle X_a,X_b\rangle = \left\{
       \begin{array}{lr}
        \frac{1}{3}\sum\limits_{i=1}^n\left(x^{-2}_{ai}+x^{-}_{ai} x^{+}_{ai}+x^{+2}_{ai}\right), & \text{if } X_a=X_b\\
       \frac{1}{4}\sum\limits_{i=1}^n\left(x^{-}_{ai}x^{-}_{bi} +x^{-}_{ai}x^{+}_{bi} +x^{+}_{ai}x^{-}_{bi}+x^{+}_{ai}x^{+}_{bi}\right), & \text{if } X_a\neq X_b
  \end{array}\right. $}
\end{equation}
From~\eqref{CIM1}, the coefficients $\beta$ can be estimated by the inverse of matrix $M$ as defined in~\eqref{CIM11}.
\begin{equation}\label{CIM11}
    \hat{\beta} = (M)^{-1} b
\end{equation}
Thus, using $\hat{\beta}$ and $\hat{\tau}_j$, the lower and upper bounds of the regressand, $Y$ can be estimated by~\eqref{CIM_lower_upper}.

The CIM method satisfies mathematical coherence of the estimated bounds of the regressand using $\tau_j$ where its value (1 or 0) depends on the sign (+ or -) of $\hat{\beta}_j$~\cite{wang2012linear}.

\subsection{The Linear Model (2015)}\label{LM_subsection}
Sun and Ralescu~\cite{sun2015linear} developed the Linear Model (LM) in 2015 based on the affine operator in the cone $\mathcal{C}=\{(x,y)\in \mathbb{R}^2|x\leq y\}$, aiming to maximise model flexibility while preserving its interpretability. This model considers both lower and upper bounds of regressors and their ranges for estimating the bounds of the regressand. In this model, $\overline{y}_i$ is expressed as a linear transformation of $\overline{x}_{ji}$ such that,
\begin{align}\centering
\begin{split}\label{LM}
y_i^-&=\sum_{j=1}^p\left(\alpha_j x_{ji}^-+\beta_j x_{ji}^+\right)+\eta+\epsilon_i^-,\\
y_i^+&=\sum_{j=1}^p\left((\alpha_j-\gamma_j) x_{ji}^-+(\beta_j+\gamma_j) x_{ji}^+\right)+\eta+\theta+\epsilon_i^+,
\end{split}
\end{align}
where the regressor coefficients, $\alpha_j,\beta_j, \eta\in \mathbb{R}$ and $\gamma_j,\theta\ge0$, $j=1,2,...,p$. Equation~\eqref{LM} can also be written in the following matrix form as $Y=X\beta+\epsilon$, where

$Y=\begin{pmatrix}
        Y^-\\
        Y^+
        \end{pmatrix}$ with $Y^-=\begin{pmatrix}
        y_1^-\\
        \vdots \\
        y_n^-
        \end{pmatrix}$ and $Y^+=\begin{pmatrix}
        y_1^+\\
        \vdots \\
        y_n^+
        \end{pmatrix}$,

\resizebox{.45\columnwidth}{!}{$X=\begin{pmatrix}
        X^* & 0\\
        X^* & X^w
        \end{pmatrix} \text{ with } X^*=\begin{pmatrix}
        1 & x_{11}^-& x_{11}^+&\dots & x_{p1}^+\\
        1 & x_{12}^-& x_{12}^+&\dots & x_{p2}^+\\
        \vdots & \vdots & \vdots&\vdots &\vdots \\
        1 & x_{1n}^-&x_{1n}^+&\dots & x_{pn}^+
        \end{pmatrix}$}
        \text{ and } $X^w=\begin{pmatrix}
        1 & x_{11}^w& \dots & x_{p1}^w\\
        1 & x_{12}^w& \dots & x_{p2}^w\\
        \vdots & \vdots &\vdots &\vdots \\
        1 & x_{1n}^w&\dots & x_{pn}^w
        \end{pmatrix}$, and  $\beta=\begin{pmatrix}
        \beta^* \\
        \beta^w
        \end{pmatrix}$ with

        $\beta^*=\begin{pmatrix}
        \eta \\
        \alpha_1\\
        \beta_1\\
        \vdots \\
        \alpha_p\\
        \beta_p
        \end{pmatrix}$ and $\beta^w=\begin{pmatrix}
        \theta \\
        \gamma_1\\
        \vdots \\
        \gamma_p
        \end{pmatrix}$.

Again, the LS estimate of the coefficients, $\hat{\beta}$ are computed by~\eqref{eq:Beta_Linear},
\begin{equation}\label{eq:Beta_Linear}
    \hat{\beta}=(X^T X)^{-1}X^TY
\end{equation}
and~\eqref{Estimated_Linear} estimates the bounds of the regressand, $Y$.
\begin{align}\label{Estimated_Linear}\centering
\begin{split}
\hat{y}_i^-&=\sum_{j=1}^p\left(\hat{\alpha}_j x_{ji}^-+\hat{\beta}_j x_{ji}^+\right)+\hat{\eta},\\
\hat{y}_i^+&=\sum_{j=1}^p\left((\hat{\alpha}_j-\hat{\gamma}_j) x_{ji}^-+(\hat{\beta}_j+\hat{\gamma}_j) x_{ji}^+\right)+\hat{\eta}+\hat{\theta}
\end{split}
\end{align}

Even though the authors assume positive constraints on range coefficients (i.e., $\gamma_j\ge0$ and $\theta\ge0$) in~\eqref{LM} to maintain mathematical coherence, the actual model does not ensure compliance with these constraints. As a result, \eqref{eq:Beta_Linear} can result in negative range coefficients ($\hat{\gamma}_j,\hat{\theta}<0$), which may lead to flipped interval bounds ($\hat{y}_i^->\hat{y}_i^+$) through \eqref{Estimated_Linear}. The authors do not discuss how to maintain these constraints in practice, though they expect that if any of the estimates of $\hat{\gamma_j}$ and $\hat{\theta}$ turns out to be negative, forcing it to be positive may lead to poor fitness of the LM model.

\subsection{The Parametrized Model (2017)}
Souza et al. \cite{souza2017parametrized} proposed the Parametrized Model (PM) in 2017, which also builds two different models for the regressand bounds. This approach extracts the best reference points from the regressors and uses them to build linear regression models for both lower and upper bounds of the regressand, unlike earlier approaches that use specific reference points on the regressors, such as center, range, interval bounds.

In this method, an interval is considered as a line segment, useful to reach all points within it. For instance, given an interval $\overline{a}$, any point $q\in\overline{a}$ can be  computed as $q=a^- (1-\lambda) + a^+\lambda$, $0\leq \lambda \leq 1$. By setting $\lambda$, $\overline{a}$ is turned into a single point. Hence, when $\lambda= 0$,  $q=a^-$ (lower bound of $\overline{a}$) and when $\lambda= 1$, $q=a^+$ (upper bound of $\overline{a}$). Similarly, $q=a^c$ (center of $\overline{a}$) when $\lambda = 0.5$. Utilizing this concept, the PM method specifies the linear regression models for the lower and upper bounds of $Y$ in~\eqref{PM_lower_upper}.
\begin{equation}\label{PM_lower_upper}
\begin{split}
    y_i^-=\beta_0^-+\sum_{j=1}^p{{\beta_j^-(1-\lambda_j) x_{ji}^-+\beta_j^-\lambda_j x_{ji}^+}}+\epsilon_i^-,\\
    y_i^+=\beta_0^++\sum_{j=1}^p{{\beta_j^+(1-\lambda_j) x_{ji}^-+\beta_j^+\lambda_j x_{ji}^+}}+\epsilon_i^+.
\end{split}
\end{equation}

Equation~\eqref{PM_lower_upper1} simplifies \eqref{PM_lower_upper} by replacing $\beta_j^-(1-\lambda_j)$ by $\alpha_j^-$ and $\beta_j^-\lambda_j$ by $\omega_j^-$ for lower bounds, and $\beta_j^+(1-\lambda_j)$ and $\beta_j^+\lambda_j$  by $\alpha_j^+$ and $\omega_j^+$  respectively.
\begin{equation}\label{PM_lower_upper1}
\begin{split}
    y_i^-=\beta_0^-+\sum_{j=1}^p{\alpha_j^- x_{ji}^-+\omega_j^- x_{ji}^+}+\epsilon_i^-\\
    y_i^+=\beta_0^++\sum_{j=1}^p{\alpha_j^+ x_{ji}^-+\omega_j^+ x_{ji}^+}+\epsilon_i^+.
\end{split}
\end{equation}

In matrix notation, the lower bound model can be expressed for all $n$ observations as $Y^- = X^*\beta^- + \epsilon^{-}$, where

\resizebox{.94\columnwidth}{!}{
  \begin{minipage}{0.9\columnwidth}
\begin{align*}\label{PM_form}
    \begin{split}
     Y^-=(y_1^-\quad y_2^-\quad\dots\quad y_n^-)^T,\quad \beta^-=(\beta_0^-\quad \alpha_1^-\quad\omega_1^-\quad\dots\quad\alpha_p^-\quad\omega_p^-)^T,\\ \epsilon^-=(\epsilon_1^-\quad \epsilon_2^-\quad\dots\quad\epsilon_n^-)^T,
     \text{ and } X^*=\begin{pmatrix}
                     1 & x_{11}^- & x_{11}^+ &\dots & x_{p1}^+\\
                     1 & x_{12}^- & x_{12}^+ &\dots & x_{p2}^+\\
                     \vdots & \vdots &\vdots &\dots &\vdots\\
                     1 & x_{1n}^- &x_{1n}^+ &\dots & x_{pn}^+\\
                     \end{pmatrix}.
\end{split}
\end{align*}
\end{minipage}}

The LS estimate of the coefficients for the lower bound model, $\beta^-$ is computed by~\eqref{eq:Beta_PM_lower}.
\begin{equation}\label{eq:Beta_PM_lower}
    \hat{\beta}^-=((X^*)^T X^*)^{-1}(X^*)^TY^-
\end{equation}

The matrix expression follows the same pattern for the upper bound model, $Y^+= X^*\beta^+ + \epsilon^{+}$, and the LS estimate of the coefficients for the upper bound model, $\beta^+$ in defined in~\eqref{eq:Beta_PM_upper}.
\begin{equation}\label{eq:Beta_PM_upper}
    \hat{\beta}^+=((X^*)^T X^*)^{-1}(X^*)^TY^+
\end{equation}
Finally, using $\hat{\beta}^-$ and $\hat{\beta}^+$, the lower and upper bounds of $Y$ are estimated using \eqref{PM_lower_upper2}.
\begin{equation}\label{PM_lower_upper2}
\begin{split}
    \hat{y}_i^-=\hat{\beta}_0^-+\sum_{j=1}^p{\hat{\alpha}_j^- x_{ji}^-+\hat{\omega}_j^- x_{ji}^+}\\
    \hat{y}_i^+=\hat{\beta}_0^++\sum_{j=1}^p{\hat{\alpha}_j^+ x_{ji}^-+\hat{\omega}_j^+ x_{ji}^+}.
\end{split}
\end{equation}

The PM method does not automatically guarantee the mathematical coherence of the bounds~\cite{souza2017parametrized}. To avoid flipping the interval bounds, the approach estimates the range of $Y$ using \eqref{eq:PM_range_estimate} before performing the regression.
\begin{equation}\label{eq:PM_range_estimate}
    \hat{Y}^w=X^*((X^*)^T X^*)^{-1}(X^*)^TY^w
\end{equation}
If all estimated ranges are positive ($\hat{y}^w\in\hat{Y}^w$), the model automatically ensures mathematical coherence. However, if at least one of the estimated ranges is negative, it applies the Box-Cox transformation~\cite{box1964analysis}, extended to intervals by the authors \cite{souza2017parametrized}, to transform the regressand so that the desirable coherence is  achieved by the PM method.

\section{Data \& Methodology}
\label{sec:methodology}
In Section II, we have discussed well-known linear regression models for intervals considering their vector/matrix representation. Table~\ref{tab:regression_models} provides a summary of these models and their features. In this section, we initially put forward an extension to the LM method \cite{sun2015linear}  which ensures mathematical coherence, thus making it more suitable for real-world data. Later, we introduce both known synthetic and real-world data sets to support the empirical exploration of the regression approaches. Lastly, a number of evaluation metrics are described through which their performance is assessed.
\begin{table*}
\centering
\caption{Summary of well-known linear regression models for interval-valued data}
\resizebox{0.97\textwidth}{!}{
\begin{tabular}[t]{llcl}
\hline
 \small{Methods} &  \small{Interval reference point(s) used} &  \small{Also applicable to numeric data} &\small{Ensure mathematical coherence on } \\
  & {\small{in the regression process}} & {\small{(interval range = 0)}}   & {\small{interval bounds}}\\\hline
CM \cite{billard2000regression}&Center&Yes&No\\
MinMax \cite{billard2002symbolic}&Lower and upper bounds&Yes&No\\
CRM \cite{neto2008centre}&Center and radius (half of range)     &No&No\\
CCRM \cite{neto2010constrained}& Center and radius (half of range)&No&Yes with positivity restrictions on range coefficients\\
CIM \cite{wang2012linear}& All points within intervals &Yes&Yes with Moore’s  linear  combination algorithm \cite{moore1966interval}\\
LM \cite{sun2015linear}  & Lower and upper bounds, and range &No&Yes with positivity restrictions as proposed in this paper\\
PM \cite{souza2017parametrized}& Any point within intervals &Yes&Yes with the interval Box-Cox transformation \cite{box1964analysis}\\\hline
 \end{tabular}}
\label{tab:regression_models}
\end{table*}

\subsection{Extension of the Linear Model [20]}
As pointed out in Section~\ref{sec:background}, the LM method~\cite{sun2015linear} considers positivity restrictions on the range coefficients to consistently estimate the regressand lower and upper bounds. However, it neither guarantees these restrictions through its model setting nor provides algorithmic solutions to maintaining them. As it is essential in practice to ensure mathematical coherence, we extend this method in the following two directions maintaining coherence.

\begin{algorithm}
\caption{A constrained LM method ($LM_c$)}
\noindent\textbf{Input:} Let $Y=\{\overline{y}_1,\overline{y}_2,...,\overline{y}_n\}$ be the regressand variable's samples and $\{X_1, X_2,..., X_p\}$ are the $p$ ($\geq 1$) be the regressor variables. \vspace{0.2cm}\\
\noindent\textbf{Output:} Let $\hat{Y}=\{\hat{\overline{y}}_1,\hat{\overline{y}}_2,...,\hat{\overline{y}}_n\}$ be the estimated $Y$
\label{algorithm_1}
\begin{algorithmic}[1]
 \State Apply the LM method at \eqref{eq:Beta_Linear} on $Y$ and $X_j$ ($1\leq j\leq p$) to estimate $\hat{\beta}^*$ and $\hat{\beta}^w$
 \If {$\hat{\theta}<0$ or $\hat{\gamma}_j<0$ where $\hat{\theta},\hat{\gamma}_j\in \hat{\beta}^w$,}
 \State Apply LHA \cite{lawson1974} to make $\hat{\theta}$ and $\hat{\gamma}_j$ positive
 \EndIf
 \State Estimate $\hat{Y}^-$ and $\hat{Y}^+$ using \eqref{Estimated_Linear}
 \State \Return $\hat{Y}=[\hat{Y}^-,\hat{Y}^+]$
\end{algorithmic}%
\end{algorithm}

\emph{(1)} We estimate a constrained LM method ($LM_c$) following the notion of the CCRM method \cite{neto2010constrained} where Lawson and Hanson’s algorithm (LHA) \cite{lawson1974} is applied to universally impose
the positivity restriction on the range coefficients, $\gamma_j$ and $\theta$. Algorithm~\ref{algorithm_1} presents this constrained LM variant, $LM_c$. Note that universally (or blindly) applying the restrictions risks sacrificing performance in terms of model-data fit in favour of ensuring mathematical coherence.

\emph{(2)} To minimise the impact on performance, we also develop a more weakly constrained LM method (LM$_w$), in line with the concept of the PM approach~\cite{souza2017parametrized} to avoid unnecessary bias in the estimation results due to universally restricting parameters. Here, we first regress the range of the regressand variable ($y_i^w$) on the range of regressors ($x_{ji}^w$) before performing the estimation using~\eqref{eq:PM_range_estimate} as the estimated range values of the regressand guarantee mathematical coherence of the estimated bounds \emph{iff} they are positive~\cite{souza2017parametrized}. Therefore, if all estimated range values are positive ($\hat{y}_i^w\geq 0, \forall i$), we compute the least squares estimation of~\eqref{eq:Beta_Linear} without any parameter restrictions. In contrast, if any of the estimated range values is negative ($\hat{y}_i^w<0, \exists i$), we put positivity restrictions on $\hat{\gamma}_j$ and $\hat{\theta}$ using the Lawson and Hanson’s algorithm (LHA)~\cite{lawson1974,neto2010constrained}. Algorithm~\ref{algorithm_2} summarises this adapted approach, LM$_w$.

\begin{algorithm}
\caption{A weakly-constrained LM method ($LM_w$)}
\noindent\textbf{Input:} Let $Y=\{\overline{y}_1,\overline{y}_2,...,\overline{y}_n\}$ be the regressand variable's samples and $\{X_1, X_2,..., X_p\}$ are the $p$ ($\geq 1$) regressor variables. \\\vspace{0.2cm}
\noindent\textbf{Output:} Let $\hat{Y}=\{\hat{\overline{y}}_1,\hat{\overline{y}}_2,...,\hat{\overline{y}}_n\}$ be the estimated $Y$
\label{algorithm_2}
\begin{algorithmic}[1]
\State Estimate the range values of regressand variable, $\hat{Y}^w=\{\hat{y}^w_1,...,\hat{y}^w_n\}$ using \eqref{eq:PM_range_estimate}
\If {all $\hat{y}^w_i\geq 0$}
   \State Apply the LM method at \eqref{eq:Beta_Linear} on $Y$ and $X_j$ ($1\leq j\leq p$) to generate $\hat{\beta}^*$ and $\hat{\beta}^w$
   \State Estimate $\hat{Y}^-$ and $\hat{Y}^+$ using \eqref{Estimated_Linear}
\Else
   \State Apply \textbf{Algorithm 1} to estimate $\hat{Y}^-$ and $\hat{Y}^+$
\EndIf
\State \Return $\hat{Y}=[\hat{Y}^-,\hat{Y}^+]$
\end{algorithmic
\end{algorithm}

\subsection{Datasets}
\label{sec:methodology.datasets}
We use both synthetic and real-world interval-valued data sets to investigate the estimation performance of the linear regression methods. The synthetic datasets are used to explore the sensitivity of the regression models in respect to different properties of the interval data,  whereas the real-world data sets are used to investigate their suitability to real-world scenarios. All data sets are described below.

\subsubsection{Synthetic Data Sets}
\label{sec:methodology.syn.datasets}
We consider a series of synthetic interval-valued datasets to investigate different aspects of model estimation. First, Set-1 and Set-2, already shown at the beginning of the paper (Fig.~\ref{fig1}), are used to examine the impact of key features (i.e., center and range). Second, an additional nine sets of intervals (Set-3 to Set-11), exhibiting different properties, are used to explore relative performance of the regression models. Specifically, we use three sets of \emph{skinny} intervals (Set-3 to Set-5) with smaller range\footnote{The `skinny' and `puffy' designs have been inspired by~Ferson et al.~\cite{ferson2007experimental}} and three \emph{puffy} sets with wider range (Set-6 to Set-8) having features of disjointedness, some overlap and full overlap. We also consider three mixed cases (Set-9 to Set-11) containing puffy and skinny intervals in separated, partially overlapping and nested states. All of these datasets are produced within the range of 0 to 10. They are given in Table~\ref{tab:dataset2}~(a)~to~(k) in Appendix including the mean ($\mu$) and standard deviation ($\sigma$) of the center and range of intervals. Figure~\ref{fig:synthetic_dataset_A} graphically presents the interval data sets, Set-3 to Set-11.

\subsubsection{Real-world Data Sets}
\label{sec:methodology.real.datasets}
We consider two different real-world interval-valued data sets, where in the first case intervals are more skinny and condensed than the second case. The first set represents medical observations; it includes interval-valued systolic and diastolic blood pressure (BP) data of 59 patients of the Hospital Valle del Nalón in Asturias, Spain~\cite{blanco2011estimation,garcia2020multiple}.
The second set considers assessment by cyber-security experts of the vulnerability of system components in a real-world scenario. Specifically, 38 cyber-security experts at the UK CESG (Communications-Electronics Security Group) assessed a range of components (referred to as hops) that are commonly encountered during a cyber-attack~\cite{miller2016modelling,ellerby2019exploring}. They rated the hops on overall difficulty for an attacker to either attack or evade them, as well as rating each on several attributes (seven for attack, and three for evade) that might affect this difficulty. The experts provided their interval-valued ratings on a scale from 0 to 100. For more details, see~\cite{miller2016modelling,ellerby2019exploring}.

\begin{figure}[b!]
\begin{subfigure}{0.23\textwidth}
\centering
\includegraphics[width=1\textwidth]{figure2a.png}
\caption{\begin{footnotesize}Set-3 (All disjoint SIs)
\end{footnotesize}}
\end{subfigure}\quad
\begin{subfigure}{0.23\textwidth}
\centering
\includegraphics[width=1\textwidth]{figure2b.png}
\caption{\begin{footnotesize}Set-4 (Some overlapping SIs)\end{footnotesize}}
\end{subfigure}\\
\begin{subfigure}{0.23\textwidth}
\centering
\includegraphics[width=1\textwidth]{figure2c.png}
\caption{\begin{footnotesize}Set-5 (All overlapping SIs)\end{footnotesize}}
\end{subfigure}\quad
\begin{subfigure}{0.23\textwidth}
\centering
\includegraphics[width=1\textwidth]{figure2d.png}
\caption{\begin{footnotesize}Set-6 (All disjoint PIs)\end{footnotesize}}
\end{subfigure}\\
\begin{subfigure}{0.23\textwidth}
\centering
\includegraphics[width=1\textwidth]{figure2e.png}
\caption{\begin{footnotesize}Set-7 (Some overlapping PIs)\end{footnotesize}}
\end{subfigure}\quad
\begin{subfigure}{0.23\textwidth}
\centering
\includegraphics[width=1\textwidth]{figure2f.png}
\caption{\begin{footnotesize}Set-8 (All overlapping PIs)\end{footnotesize}}
\end{subfigure}\\
\begin{subfigure}{0.23\textwidth}
\centering
\includegraphics[width=1\textwidth]{figure2g.png}
\caption{\begin{footnotesize}Set-9 (Disjoint MIs)\end{footnotesize}}
\end{subfigure}\quad
\begin{subfigure}{0.23\textwidth}
\centering
\includegraphics[width=1\textwidth]{figure2h.png}
\caption{\begin{footnotesize}Set-10 (Some overlapping MIs)\end{footnotesize}}
\end{subfigure}\\
\begin{subfigure}{0.46\textwidth}
\centering
\includegraphics[width=0.5\textwidth]{figure2i.png}
\caption{\begin{footnotesize}Set-11 (Nested/all overlapping MIs)\end{footnotesize}}
\end{subfigure}
\caption{2D visualization of synthetic datasets, Set-3 to Set-11. SI=skinny interval, PI=puffy interval, and MI=mixed interval.}
\label{fig:synthetic_dataset_A}
\end{figure}

\subsection{Evaluation Metrics}
\label{sec:methodology.metrics}
We use the following three metrics to evaluate the performance of the linear regression models. Each of the evaluation metrics is defined in respect to the regressand (samples), $Y=\{\overline{y}_1,\overline{y}_2,...,\overline{y}_n\}$ and its estimated value(s),  $\hat{Y}=\{\hat{\overline{y}}_1,\hat{\overline{y}}_2,...,\hat{\overline{y}}_n\}$.

\noindent\emph{1) Root Mean Squared Error (RMSE)} estimates error by considering the average deviation of the estimated values from the actual ones~\cite{neto2010constrained}. The \emph{RMSE} for the lower and upper bounds of $Y$ is defined as,
\begin{equation*}\label{eq:RMSE}
\resizebox{0.94\columnwidth}{!}{$
    RMSE^-=\sqrt{\frac{1}{n}\sum\limits_{i=1}^{n}(y_i^--\hat{y}_i^-)^2} \text{ and }
    RMSE^+=\sqrt{\frac{1}{n}\sum\limits_{i=1}^{n}(y_i^+-\hat{y}_i^+)^2},
$}
\end{equation*}
where
a smaller \emph{RMSE} value indicates a better fitted model.

\vspace{0.1mm}\noindent\emph{2) Mean Absolute Error (MAE)} calculates the average of absolute difference between the actual and estimated values~\cite{iccen2016error}. For the upper and lower bounds of $Y$, the \emph{MAE} is defined as,

\begin{equation*}\label{eq:MAE}
\resizebox{0.94\columnwidth}{!}{ $
MAE^-=\frac{1}{n}\sum\limits_{i=1}^n\left|y_i^--\hat{y}_i^-\right| \text{ and } MAE^+=\frac{1}{n}\sum\limits_{i=1}^n\left|y_i^+-\hat{y}_i^+\right|,$}
\end{equation*}
where $0\leq MAE^-,MAE^+\leq\mathbb{R}^+$. Again, a lower value of \emph{MAE} implies a better fitted model.

\vspace{0.1mm}\noindent\emph{3) Mean Magnitude of Relative Error (MMRE)} measures the mean of dispersion of the actual and estimated values in respect to the actual values~\cite{fagundes2016quantile}. It is calculated as follows:
\begin{equation*} \label{eq:MMRE}
    MMRE=\frac{1}{n}\sum\limits_{i=1}^n\frac{1}{2}\left\{\left|\frac{y^-_i-\hat{y}^-_i}{y^-_i}\right|+\left|\frac{y^+_i-\hat{y}^+_i}{y^+_i}\right|\right\},
\end{equation*}
where $0\leq MMRE\leq \mathbb{R}^+$. In this case, the lower the \emph{MMRE}, the better the model.

\section{Results}
\label{sec:results}
We explore the performance of the existing regression models together with the two proposed variants of the LM method (LM$_c$ and LM$_w$) for both synthetic and real-world interval data sets (see Section~\ref{sec:methodology.datasets}). Their performance is assessed through different metrics as described in Section~\ref{sec:methodology.metrics}.

\subsection{Performance Analysis with Synthetic Data Sets}
\label{subsec:syn_results}
In this section, we explore following two aspects to investigate the sensitivity of different regression models in respect to various features of intervals. They are as follows:
\begin{itemize}[leftmargin=*]
\item Impact of key features (center and range) of the regressor on those of the regressand.
\item Performance of regression methods in respect to different types of interval-valued datasets.
\end{itemize}

\subsubsection {Impact of key features (center and range) of regressor on those of regressand}
Here, we consider Set-1 and Set-2 (see Fig.~1~(a)~and~(b)) which we introduced at the beginning of the paper. As noted, Set-1 shows a commonly considered and intuitive scenario where uncertainty and position in regressor and regressand vary in unison. On the other hand, Set-2 captures an extreme relationship with the range of the regressand independent of the range of the regressor. Incidentally, it is straightforward to consider the mirrored (`horizontal') extreme case, but we do not consider it here for space limitations.

We apply the regression models for both Sets 1 and 2 to estimate the coefficients (shown in Appendix). Tables~\ref{tab:my_label1}~and~\ref{tab:my_label2} present their fitness. Beyond considering model performance, in order to clearly communicate and visualise the behaviour of the regression models, we introduce 3D surface plots which capture the (degree of) change in the regressand's key features (center, range) for given changes in the regressor's key features (center, range). For conciseness, we present combined (for center and range) 3D surface plots for all methods in Appendix (see Fig.~\ref{fig:case-1}~and~\ref{fig:case-2}). Here, for clarity and illustration, we independently present the impact of the center and range of the test regressor on the center (subfigures (a)) and the range (subfigures (b)) of the regressand, for the CRM and PM methods in Fig.~\ref{fig:synthetic2-case1} and~\ref{fig:synthetic2-case11} and Fig.~\ref{fig:synthetic2-case2} and~\ref{fig:synthetic2-case22}, respectively.

Figure~\ref{fig:synthetic2-case1}~(a)~and~(b) show that for the CRM method, the estimated center of Set-1's regressand is influenced only by the center of the regressor, and the estimated range of Set-1 is influenced only by the range of the regressor. Thus, the CRM method is not in line with our intuition considering the dataset. On the other hand, Fig.~\ref{fig:synthetic2-case11}~(a)~and~(b) show that both changes in the estimated center and range of the regressand of Set-1 using the PM method are associated with changes in the center and range of the regressor. Thus, the PM method performs according to our expectation. Further, the MinMax method is found to be performing in-line with the PM method. However, the rest of the methods, like the CRM, do not (see Appendix, Fig.~\ref{fig:case-1}).

\begin{figure}[t!]
\begin{subfigure}{0.24\textwidth}
\includegraphics[width=0.985\textwidth]{figure3a.png}
\caption{In respect to $\hat{Y}^c$}
\end{subfigure}
\begin{subfigure}{0.24\textwidth}
\includegraphics[width=0.985\textwidth]{figure3b.png}
\caption{In respect to $\hat{Y}^w$ ($\simeq 2\hat{Y}^r$)}
\end{subfigure}
\caption{Regression behaviour of the CRM method for Set-1.
\label{fig:synthetic2-case1}
\end{figure}
\begin{figure}[t!]
\begin{subfigure}{0.24\textwidth}
\includegraphics[width=0.985\textwidth]{figure4a.png}
\caption{In respect to $\hat{Y}^c$}
\end{subfigure}
\begin{subfigure}{0.24\textwidth}
\includegraphics[width=0.985\textwidth]{figure4b.png}
\caption{In respect to $\hat{Y}^w$ ($\simeq 2\hat{Y}^r$)}
\end{subfigure}
\caption{Regression behaviour of the PM method for Set-1.
\label{fig:synthetic2-case11}
\end{figure}
\begin{table}[t!]
    \centering
\caption{Regression performance for Set-1}
\resizebox{0.9\hsize}{!}{
\begin{tabular}{cccccc} \hline
Methods & $RMSE^-$ & $RMSE^+$ & $MAE^-$ & $MAE^+$ & $MMRE$
\\\hline
CM &	1.505&	1.51&	1.354&	1.354&	0.427\\
MinMax&	0.002&	0.113&	0.001&	0.109&	0.005\\
CRM	&0.16&	0.055&	0.14&	0.047&	0.0333\\
CCRM&	0.251&	0.227&	0.218&	0.198&	0.055\\
CIM	&1.503&	1.517&	1.356&	1.356&	0.427\\
PM	&\textbf{0.0}&	\textbf{0.04}&	\textbf{0.0}&	\textbf{0.031}&	\textbf{0.001}\\
LM$_c$	&0.097&	0.428&	0.085&	0.381&	0.034\\
LM$_w$	&0.097&	0.105&	0.085&	0.088&	0.024\\\hline
\begin{scriptsize}\textbf{Bold} = Best\end{scriptsize}&&&&&
\end{tabular}}
\label{tab:my_label1}
\end{table}

Figure~\ref{fig:synthetic2-case2}~(a) visualises how the estimated center of Set-2's regressand, using the CRM method, is positively influenced by changes in the center, but not the range, of the regressor -- as expected. However, Figure~\ref{fig:synthetic2-case2}~(b) illustrates how the range of Set-2 remains mainly unexplained by characteristics of the regressor using the CRM method.
By contrast, Fig.~\ref{fig:synthetic2-case22}~(a)~and~(b) show how when using the PM method, both estimated center and range of the regressand of Set-2 are well-explained by the center of the regressor -- in line with our expectations. From Fig.~\ref{fig:case-2} (Appendix), we can observe that the MinMax method shows similar performance to the PM method. However, the remaining methods again do not meet intuitive expectations.

\subsubsection {Broader performance of regression methods}
Here we consider nine synthetic cases of Fig.~\ref{fig:synthetic_dataset_A}
(discussed in Section~\ref{sec:methodology.syn.datasets}). Table~\ref{tab:synthetic_dataset_A_Table}~(a)~--~(i) presents the fitness of different regression methods for each of these datasets.
We summarize the results
below, before discussing performance more broadly -- across all datasets (Set-3 to Set-11).

\begin{figure}[t!]
\begin{subfigure}{0.24\textwidth}
\includegraphics[width=0.99\textwidth]{figure5a.png}
\caption{In respect to $\hat{Y}^c$}
\end{subfigure}
\begin{subfigure}{0.24\textwidth}
\includegraphics[width=0.99\textwidth]{figure5b.png}
\caption{In respect to $\hat{Y}^w$ ($\simeq 2\hat{Y}^r$)}
\end{subfigure}
\caption{Regression behaviour of the CRM method for Set-2.}
\label{fig:synthetic2-case2}
\end{figure}
\begin{figure}[t!]
\begin{subfigure}{0.24\textwidth}
\includegraphics[width=0.99\textwidth]{figure6a.png}
\caption{In respect to $\hat{Y}^c$}
\end{subfigure}
\begin{subfigure}{0.24\textwidth}
\includegraphics[width=0.99\textwidth]{figure6b.png}
\caption{In respect to $\hat{Y}^w$ ($\simeq 2\hat{Y}^r$)}
\end{subfigure}
\caption{Regression behaviour of the PM method for Set-2.}
\label{fig:synthetic2-case22}
\end{figure}
\begin{table}[t!]
\centering
\caption{Regression performance for Set-2}
\resizebox{0.9\hsize}{!}{
\begin{tabular}{cccccc} \hline
Methods & $RMSE^-$ & $RMSE^+$ & $MAE^-$ & $MAE^+$ & $MMRE$ \\\hline
CM	&1.65	&1.657	&1.491	&1.491	&0.467\\
MinMax	&\textbf{0.002}	&\textbf{0.153}	&\textbf{0.001}	&0.148	&\textbf{0.006}	\\
CRM	&0.575	&0.592	&0.53	&0.535	&0.138	\\
CCRM	&0.703	&0.719	&0.61	&0.625	&0.167	\\
CIM	&1.65	&1.657	&1.491	&1.491	&0.467\\
PM	&\textbf{0.002}	&\textbf{0.153}	&\textbf{0.001}	&\textbf{0.145}	&\textbf{0.006}	\\
LM$_c$	&0.573	&1.007	&0.528	&0.929	&0.156\\
LM$_w$	&0.573	&0.593	&0.528	&0.537	&0.138\\\hline
\begin{scriptsize}\textbf{Bold} = Best\end{scriptsize}&&&&&
\end{tabular}}
\label{tab:my_label2}
\end{table}

\begin{table*
\caption{Performance of different regression models for interval-valued synthetic datasets (Set-3 to Set-11)}
\begin{subtable}{.33\linewidth}
\centering
\caption{Set-3 (All disjoint skinny intervals)}
\renewcommand{\arraystretch}{1.1}
\resizebox{0.95\columnwidth}{!}{
\begin{Huge}
\begin{tabular}{ccccccc} \hline
Methods & $RMSE^-$ & $RMSE^+$ & $MAE^-$ & $MAE^+$ & $MMRE$ \\\hline
CM&1.215	&1.181	&0.972	&1.035	&0.04\\
MinMax&1.204	&1.169	&1.004	&1.002	&0.04\\
CRM&1.185	&1.189	&0.986	&1.021	&0.04\\
CCRM&1.195	&1.179	&0.995	&1.011	&0.04\\
CIM&1.215	&1.181	&0.972	&1.035	&0.04\\
PM&\textbf{0.343}	&\textbf{0.355}	&\textbf{0.291}	&0.312	&\textbf{0.016}\\
LM$_c$&\textbf{0.343}	&0.356	&0.293	&0.317	&\textbf{0.016}\\
LM$_w$&\textbf{0.343}	&\textbf{0.355}	&0.293	&\textbf{0.310}	&\textbf{0.016}\\\hline\vspace{0.1mm}
\end{tabular}\end{Huge}}
\end{subtable}%
\begin{subtable}{.33\linewidth}
\centering
\caption{Set-4  (Some overlapping  skinny  intervals)}
\renewcommand{\arraystretch}{1.1}
\resizebox{0.95\columnwidth}{!}{
\begin{Huge}
\begin{tabular}{ccccccc} \hline
Methods & $RMSE^-$ & $RMSE^+$ & $MAE^-$ & $MAE^+$ & $MMRE$ \\\hline
CM&0.656	&0.601	&0.592	&0.407	& 0.021	\\
MinMax&0.631&	0.574&	0.483&	0.432&	0.015	\\
CRM&0.616&	0.591&	0.475&	0.451&	0.016\\
CCRM&0.626&	0.581&	0.482&	0.443&	0.016\\
CIM&0.656&	0.601	&0.592&	0.407&	0.021\\
PM&\textbf{0.263}&\textbf{0.192}&	\textbf{0.249}&	\textbf{0.181}&	\textbf{0.014}\\
LM$_c$&\textbf{0.263}&	0.193&\textbf{0.249}&	0.185&\textbf{0.014}\\
LM$_w$&\textbf{0.263}&	\textbf{0.192}&	\textbf{0.249}&	\textbf{0.181}&	\textbf{0.014}\\\hline\vspace{0.1mm}
\end{tabular}\end{Huge}}
\end{subtable}%
\begin{subtable}{.33\linewidth}
\centering
\caption{Set-5 (All overlapping skinny intervals)}
\renewcommand{\arraystretch}{1.1}
\resizebox{0.95\columnwidth}{!}{
\begin{Huge}
\begin{tabular}{ccccccc} \hline
Methods & $RMSE^-$ & $RMSE^+$ & $MAE^-$ & $MAE^+$ & $MMRE$ \\\hline
CM&	0.370&	0.362&	0.341&	0.341&	0.021\\
MinMax&	0.141&	\textbf{0.115}&\textbf{0.120}&	\textbf{0.107}&	\textbf{0.006}\\
CRM&	0.143&	0.120&	0.125&	0.108&	\textbf{0.006}\\
CCRM&	0.143&	0.121&	0.128&	0.111&	\textbf{0.006}\\
CIM&	0.204&	0.177&	0.159&	0.149&	0.009\\
PM&	\textbf{0.139}&	\textbf{0.115}&	0.128&	\textbf{0.107}&	\textbf{0.006}\\
LM$_c$&	0.143&	0.122&	0.126&	0.113&	\textbf{0.006}\\
LM$_w$&	0.143&	0.120&	0.126&	0.109&	\textbf{0.006}\\\hline\vspace{0.1mm}
\end{tabular}\end{Huge}}
\end{subtable}\\%
\begin{subtable}{.33\linewidth}
\centering
\caption{Set-6 (All disjoint puffy intervals)}
\renewcommand{\arraystretch}{1.1}
\resizebox{0.95\columnwidth}{!}{
\begin{Huge}
\begin{tabular}{ccccccc} \hline
Methods & $RMSE^-$ & $RMSE^+$ & $MAE^-$ & $MAE^+$ & $MMRE$ \\\hline
CM&	1.996&	1.999&	1.796&	1.568&	0.140\\
MinMax&	1.681&	1.685&	1.299&	1.396&	0.107\\
CRM&	1.692&	1.673&	1.314&	1.389&	0.107\\
CCRM&	1.679&	1.687&	1.297&	1.391&	0.106\\
CIM&	1.976&	1.986&	1.777&	1.564&	0.139\\
PM&	\textbf{1.456}&	\textbf{1.344}&	\textbf{1.143}&	1.061&	\textbf{0.059}\\
LM$_c$&	1.460&	1.350&	1.147&	1.059&	0.061\\
LM$_w$&	1.460&	1.349&	1.147&	\textbf{1.057}&	\textbf{0.059}\\\hline\vspace{0.1mm}
\end{tabular}\end{Huge}}
\end{subtable}%
\begin{subtable}{.33\linewidth}
\centering
\caption{Set-7 (Some overlapping puffy intervals)}
\renewcommand{\arraystretch}{1.1}
\resizebox{0.95\columnwidth}{!}{
\begin{Huge}
\begin{tabular}{ccccccc} \hline
Methods & $RMSE^-$ & $RMSE^+$ & $MAE^-$ & $MAE^+$ & $MMRE$ \\\hline
CM&	1.751&	1.697&	1.673&	1.220&	0.094\\
MinMax&	1.434&	1.367&	1.070&	1.034&	\textbf{0.052}\\
CRM&	1.442&	1.357&	1.1067&	1.036&	0.053\\
CCRM&	1.437&	1.364&	1.090&	1.038&	0.053\\
CIM&	1.753&	1.698&	1.674&	1.220&	0.094\\
PM&	\textbf{1.411}&	\textbf{1.283}&	\textbf{1.067}&	0.984&	0.053\\
LM$_c$&	1.415&	1.289&	1.071&	0.983&	0.055\\
LM$_w$&	1.415&	1.288&	1.071&	\textbf{0.981}&	0.054\\\hline\vspace{0.1mm}
\end{tabular}\end{Huge}}
\end{subtable}%
\begin{subtable}{.33\linewidth}
\centering
\caption{Set-8 (All overlapping puffy intervals)}
\renewcommand{\arraystretch}{1.1}
\resizebox{0.95\columnwidth}{!}{
\begin{Huge}
\begin{tabular}{ccccccc} \hline
Methods & $RMSE^-$ & $RMSE^+$ & $MAE^-$ & $MAE^+$ & $MMRE$ \\\hline
CM&	1.308&	1.315&	1.238&	1.238&	0.075\\
MinMax&	0.416&	0.437&	0.387&	0.398&	0.022\\
CRM&	0.416&	0.445&	0.365&	0.387&	0.020\\
CCRM&	0.410&	0.452&	0.348&	0.404&	0.020\\
CIM&	1.073&	1.084&	0.986&	0.986&	0.059\\
PM&	\textbf{0.396}&	\textbf{0.406}&	0.321&	\textbf{0.328}&	0.017\\
LM$_c$&	0.415&	0.429&\textbf{0.313}&	0.369&	\textbf{0.016}\\
LM$_w$&	0.415&	0.425&	\textbf{0.313}&	0.348&	\textbf{0.016}\\\hline\vspace{0.1mm}
\end{tabular}\end{Huge}}
\end{subtable}\\%
\begin{subtable}{.33\linewidth}
\centering
\caption{Set-9 (Disjoint mixed intervals)}
\renewcommand{\arraystretch}{1.1}
\resizebox{0.95\columnwidth}{!}{
\begin{Huge}
\begin{tabular}{ccccccc} \hline
Methods & $RMSE^-$ & $RMSE^+$ & $MAE^-$ & $MAE^+$ & $MMRE$ \\\hline
CM&	1.871&	2.449&	1.616&	1.882&	0.207\\
MinMax&	1.596&	2.236&	1.349&	1.686&	0.152\\
CRM&	1.886&	1.911&	1.497&	1.526&	0.152\\
CCRM&	1.886&	1.911&	1.497&	1.526&	0.152\\
CIM&	1.828&	2.298&	1.550&	1.796&	0.195\\
PM&	\textbf{1.210}&\textbf{1.238}&\textbf	{0.065}&	1.089&	\textbf{0.084}\\
LM$_c$&	\textbf{1.210}&\textbf{1.238}&	1.066&	\textbf{1.087}&	0.085\\
LM$_w$&	\textbf{1.210}&\textbf{1.238}&	1.066&	\textbf{1.087}&	0.085\\\hline
\end{tabular}\end{Huge}}
\end{subtable}%
\begin{subtable}{.33\linewidth}
\centering
\caption{Set-10 (Some overlapping mixed intervals)}
\renewcommand{\arraystretch}{1.1}
\resizebox{0.95\columnwidth}{!}{
\begin{Huge}
\begin{tabular}{ccccccc} \hline
Methods & $RMSE^-$ & $RMSE^+$ & $MAE^-$ & $MAE^+$ & $MMRE$ \\\hline
CM&	1.361&	1.572&	1.293&	1.307&	0.090\\
MinMax&	0.994&	1.235&	0.906&	0.891&	0.044\\
CRM&	1.107&	1.105&	0.938&	0.918&	\textbf{0.041}\\
CCRM&	1.107&	1.105&	0.938&	0.918&	\textbf{0.041}\\
CIM	&1.344&	1.549&	1.278&	1.281&	0.087\\
PM	&\textbf{0.988}     & \textbf{0.977}&	\textbf{0.883}&	0.869&	0.055\\
LM$_c$ & \textbf{0.988} & \textbf{0.977}&	0.886&	\textbf{0.866} &	0.055\\
LM$_w$	&\textbf{0.988}	& \textbf{0.977}&   0.886&  \textbf{0.866} &0.055\\\hline
\end{tabular}\end{Huge}}
\end{subtable}%
\begin{subtable}{.33\linewidth}
\centering
\caption{Set-11 (Nested mixed intervals)}
\renewcommand{\arraystretch}{1.1}
\resizebox{0.95\columnwidth}{!}{
\begin{Huge}
\begin{tabular}{ccccccc} \hline
Methods & $RMSE^-$ & $RMSE^+$ & $MAE^-$ & $MAE^+$ & $MMRE$ \\\hline
CM&	1.008&	1.022&	0.926&	0.926&	0.104\\
MinMax&	0.048&	0.048&	0.046&	0.045&	0.004\\
CRM&	0.027&	0.028&	0.019&	0.026&	\textbf{0.002}\\
CCRM&	0.027&	0.028&	0.019&	0.026&	\textbf{0.002}\\
CIM&	0.945&	0.958&	0.870&	0.870&	0.097\\
PM&	\textbf{0.018}&	\textbf{0.020}&	0.014&	\textbf{0.016}&	\textbf{0.002}\\
LM$_c$&	\textbf{0.018}&	\textbf{0.020}&	\textbf{0.013}&	\textbf{0.016}&	\textbf{0.002}\\
LM$_w$&	\textbf{0.018}&	\textbf{0.020}&	\textbf{0.013}&	\textbf{0.016}&	\textbf{0.002}\\\hline
\end{tabular}\end{Huge}}
\end{subtable}\vspace{0.5mm}\\
\qquad\begin{scriptsize}\textbf{Bold} = Best\end{scriptsize}
\label{tab:synthetic_dataset_A_Table}
\end{table*}

\begin{itemize}[leftmargin=*]
\item \emph{Set-3 (All disjoint skinny intervals)}:
The PM and LM$_w$ methods produce the best fit for all disjoint skinny intervals. The LM$_c$ method also equally fits well in respect to two metrics ($RMSE^-$ and $MMRE$).

\item \emph{Set-4 (Some overlapping skinny intervals)}:  Again, the PM and LM$_w$ methods produce the best fit. The LM$_c$ method also gives the same fitness in respect to three metrics ($RMSE^-$, $MAE^-$, and $MMRE$). The fitness by other methods are improved with some overlapping of intervals.

\item \emph{Set-5 (All overlapping skinny intervals)}:
The PM and MinMax methods come out as the best approach. Among others, CRM, CCRM, LM$_c$ and LM$_w$ methods produce quite similar results in respect to some evaluation metrics.

\item \emph{Set-6 (All disjoint puffy intervals)}:
The PM method gives the best fit, though the LM$_c$ and LM$_w$ methods also produce quite similar results for the same feature of puffy intervals.

\item \emph{Set-7 (Some overlapping puffy intervals)}:
The PM, MinMax and LM$_w$ emerge as the best approaches in respect to different metrics, though the CRM, CCRM and LM$_c$ methods also produce quite similar good fitness as to the $MMRE$.

\item \emph{Set-8 (All overlapping puffy intervals)}: For all overlapping puffy intervals, the PM method gives the best fit in respect to three metrics and the LM$_w$ method is for two evaluation metrics. Except CM and CIM methods, growing overlap among intervals appears to lower estimation errors by the CRM, CCRM and MinMax methods.

\item \emph{Set-9 (Disjoint mixed intervals)}:
The LM$_w$, LM$_c$ and PM methods perform the best among all. All other methods give relatively low fitness. In particular, the CM and CIM methods achieve the worst results.

\item \emph{Set-10 (Some overlapping mixed intervals)}:
The PM, LM$_c$ and LM$_w$ methods bring about quite similar model fit for almost all metrics. The CRM and CCRM methods give better fitness results in respect to one metric ($MMRE$).

\item \emph{Set-11 (Nested/all overlapping mixed intervals)}:
The LM$_c$, LM$_w$, and PM methods give the best fit. The CRM and CCRM methods also closely fit the data as to $MMRE$.
\end{itemize}

Overall, the PM, and LM$_w$ (often LM$_c$ also) methods perform best and give reliable -- in fact almost identical -- fitness results for disjoint skinny or puffy intervals. Although with growing overlap among intervals (skinny or puffy), model fitness by the CRM, CCRM and MinMax methods are improved, they still do not consistently reach the performance of the PM and LM methods. For the mixed sets with puffy and skinny intervals, where they are disjoint or exhibit some overlap, both PM and two variants of LM methods give similar (best) fit. Other methods, particularly the CRM and CCRM, give nearly identical results, in respect to some metrics, for the all overlapping (nested) mixed case. In all cases, the estimates by the CIM and CM methods deviate substantially from the actual values and hence they come out as poor performers.

As well as the empirical analyses reported above, we carried out additional experiments going beyond the scope of this paper -- exploring the impact of differences in mean interval size/range, differences in the standard deviation of the mean size/range, and the variation in interval centers (i.e., dispersion), on model fitness or estimation errors. Initial results suggest that estimation errors tend to grow for all methods in respect to increasing dispersion of interval ranges (as expected). Further, greater dispersion of centers for relatively scattered interval sets increases estimation errors. We will extend  this initial investigation in more detail in the future.

\subsection{Real-world Data Sets}
\label{subsec:real_results}
In this section, the two different real-world data sets described in Section~\ref{sec:methodology.datasets} are used to evaluate the performance of the different regression methods. The fitness of these methods are again compared using the same evaluation metrics (i.e., \emph{RMSE}, \emph{MAE}, \emph{MMRE}) as before.

\subsubsection{Systolic and Diastolic Blood Pressure (BP) Data Set}
For this data set, we regress the interval-valued diastolic BP (regressand) on those of the systolic BP (regressor) to estimate the bounds of diastolic BP. Figure~\ref{fig:systolic} presents this data set and the performance of regression methods using it.

Figure~\ref{fig:systolic}~(b) shows that the PM, LM$_w$, and LM$_c$ methods produce comparable performance in estimating the regressand bounds. In respect to some evaluation metrics (\emph{RMSE$^+$} and \emph{RMSE$^-$}, and \emph{MAE$^+$}), the PM produces the best estimation, whereas for other metrics (\emph{MAE$^-$}, and \emph{MMRE}) the LM$_w$ and LM$_c$ methods come out to be the best. The CRM and CCRM methods also equally best fit the data according to the \emph{MMRE} metric. In this case, the CM and CIM methods give worse fit.

\begin{figure}
\centering
\begin{subfigure}{0.42\textwidth}
\centering
\includegraphics[width=0.95\textwidth]{figure7.png}
\caption{View of blood pressure data set \protect\cite{blanco2011estimation} }
\end{subfigure}
\begin{subfigure}{0.49\textwidth}
\vspace{1mm}
\centering
\resizebox{0.96\hsize}{!}{
\begin{tabular}{ccccccc} \hline
\large{Methods} & $RMSE^-$ & $RMSE^+$ & $MAE^-$ & $MAE^+$ & $MMRE$ \\\hline
CM	& 1.259 &	1.617 &1.097&	1.395&	0.006\\
MinMax	&0.793&	1.275&	0.651&	1.058&	0.002\\
CRM	&0.881&	1.256&	0.653&	1.029&	\textbf{0.001}\\
CCRM	&0.881	&1.256&	0.653&	1.029&	\textbf{0.001}\\
CIM		&2.104	&2.249	&1.951	&1.907	&0.011	\\
PM	&\textbf{0.780}	&\textbf{1.208}	&0.638	&\textbf{0.999}	&0.002	\\
LM$_c$	&0.784 &1.211	&\textbf{0.631}	&1.008	&\textbf{0.001}	\\
LM$_w$	&0.784	&1.211	&\textbf{0.631}	&1.008	&\textbf{0.001}	\\\hline
\begin{scriptsize}\textbf{Bold} = Best\end{scriptsize}&&&&&
\end{tabular}}
\caption{Model Performance
\end{subfigure}
\caption{Graphical presentation of systolic and diastolic blood pressure data set \protect\cite{blanco2011estimation} and performance of different regression models as to evaluation metrics using this set. }
\label{fig:systolic}
\end{figure}

\subsubsection{Cyber-security Data Set}
\label{subsec:data}
Section~\ref{sec:methodology.datasets} already provided an overview of this data set. Here, we only consider the evade hops as with three attributes they provide a more concise basis for presenting the results within the constraints of this method-centric paper. For each evade hop, 38 experts gave ratings on following three attributes: complexity (\emph{C}), availability of information (\emph{A}) and maturity (\emph{M}). They were also asked to provide an overall rating (\emph{O}) on how difficult it would be \emph{overall} for an attacker to evade a given hop. These ratings were given on a scale from 0 to 100. Table~\ref{tab:evade} summarises the attributes for the evade hops.

We explore how well the regression methods estimate the interval-valued \emph{overall} difficulty ratings for selected evade components, based upon the three interval-valued attribute ratings, and by comparison with experts' actual overall ratings. For conciseness, we present two evade hops -- number 4 and number 22, shown in Fig.~\ref{fig:real_4} and~\ref{fig:real_22} respectively.

\begin{table}
\centering
\caption{Attributes for evade hops, the questions used in the study, and the variables' name used in the analysis}
\renewcommand{\arraystretch}{1}
\begin{tabular}{c l l} \hline
Variable & Attributes & Related Question\\\hline
$C$ & Complexity & How complex is the job of providing this\\
  &            &  kind of defence?\\$A$ & Availability of & How likely is that there will be publicly \\
 & information & available information that could help with\\
 &  & evading defence?\\
$M$ & Maturity & How mature is this type of technology?\\\hline
$O$ & Overall rating & Overall, how difficult would it be for an \\
  &                & attacker to do this?\\\hline
\end{tabular}\vspace{-1.0em}
\label{tab:evade}
\end{table}

Figure~\ref{fig:real_4}~(b) shows that the PM method appears to be the best performer for Evade-4 in respect to four evaluation metrics (\emph{RMSE$^-$}, \emph{RMSE$^+$}, \emph{MAE$^+$}, and \emph{MMRE}) whereas the LM$_w$ method produces the joint best fit on two evaluation metrics (\emph{MAE$^-$} and \emph{MMRE}). The LM$_c$ method also performs very close to the LM$_w$ method. Again, the CM method provides poor fit to the data. As for the Evade-4 data set, Fig.~\ref{fig:real_22}~(b) displays similar results for Evade-22.  In particular, Fig.~\ref{fig:real_22} reveals the PM method as the best approach by most metrics. Further, the LM$_w$, LM$_c$, and CCRM methods come out as suitable measures in respect to either $MAE^+$ or $MMRE$ metric. We provide all estimated coefficients for both synthetic and real-world data sets in Appendix (see Table~\ref{tab:my_label11} and \ref{tab:my_label22}).

\begin{figure}
\centerin
\begin{subfigure}{0.42\textwidth}
\centering
\includegraphics[width=0.94\textwidth]{figure8.png}
\caption{View of Evade-4}
\end{subfigure}
\begin{subfigure}{0.49\textwidth}
\vspace{1mm}
\centering
\resizebox{0.96\hsize}{!}{
\begin{tabular}{ccccccc} \hline

\large{Methods} & $RMSE^-$ & $RMSE^+$ & $MAE^-$ & $MAE^+$ & $MMRE$ \\\hline
CM	&24.071	&25.124	&20.745	&20.478	&0.365	\\
MinMax	&21.737	&22.856	&17.948	&19.946	&0.292	\\
CRM	&22.269	&21.562	&18.631	&18.641	&0.30	\\
CCRM	&22.164	&21.756	&18.385	&18.82	&0.288	\\
CIM		&22.593	&22.582	&19.297	&19.158	&0.325		\\
PM	& \textbf{21.681}	&\textbf{20.984}	&17.852	&\textbf{17.541}	&\textbf{0.285}	\\
LM$_c$	&21.714	&21.194	&\textbf{17.807}&17.909	&0.286	\\
LM$_w$	&21.714	&21.019	&\textbf{17.807}	&17.671	&\textbf{0.285}	\\\hline
\begin{scriptsize}\textbf{Bold} = Best\end{scriptsize}&&&&&
\end{tabular}
\caption{Model Performance}
\end{subfigure}
\caption{Graphical presentation of Evade-4 data set \protect\cite{miller2016modelling} and performance of different regression models as to evaluation metrics using this set.}
\label{fig:real_4}
\end{figure}
\begin{figure}[t!]
\centering
\begin{subfigure}{0.42\textwidth}
\centering
\includegraphics[width=0.94\textwidth]{figure9.png}
\caption{View of Evade-22}
\end{subfigure}
\begin{subfigure}{0.49\textwidth}
\vspace{2mm}\centering
\resizebox{0.96\hsize}{!}{
\begin{tabular}{ccccccc} \hline
\large{Methods} & $RMSE^-$ & $RMSE^+$ & $MAE^-$ & $MAE^+$ & $MMRE$ \\\hline
CM	&27.132	&28.312	&22.759	&23.112	&0.393	\\
MinMax	&14.449	&17.972	&12.051	&14.591	&0.172	\\
CRM	&12.491	&15.399	&10.194	&12.703	&0.099	\\
CCRM	&12.487	&15.403	&10.186	&12.714	&\textbf{0.097}	\\
CIM		&15.636	&18.706	&13.405	&16.011	&0.251	\\
PM	&\textbf{12.333}	&\textbf{15.22}	&\textbf{9.862}	&12.756	& 0.109	\\
LM$_c$	&12.373	&15.255	&9.903	& \textbf{12.652}	&0.112	\\
LM$_w$	&12.373	&15.253	&9.903	&\textbf{12.652}	&0.112	\\\hline
\begin{scriptsize}\textbf{Bold} = Best\end{scriptsize}&&&&&
\end{tabular}}
\caption{Model Performance}
\end{subfigure}
\caption{Graphical presentation of Evade-22 data set \protect\cite{miller2016modelling} and performance of different regression models as to evaluation metrics using this set.}
\label{fig:real_22}
\end{figure}

\section{Discussion on Suitability of Existing Regression Approaches}
\label{sec:discussion}
This section discusses the suitability of different regression techniques, in respect to interval data sets with different properties (e.g., skinny, puffy or mixed), and based upon the results of the analyses conducted above. As highlighted, the performance of different approaches is sensitive to the characteristics and relative placement of the intervals to which they are applied.

The experiments consistently show that the PM and LM$_w$ (and often LM$_c$) are the generally best fitting models -- with reliably good performance for skinny/puffy or mixed sets of intervals, and when they are disjoint or overlapping. With increasing overlap among intervals, estimation errors decrease for the CRM, CCRM and MinMax methods, but not so far as to equal the performance of either PM or LM methods. In all cases, both MinMax and CRM require checking as to whether they maintain consistency of interval bounds. The real-world data experiments largely substantiate the synthetic experiments in terms of relative model performance.

Results reveal that the LM$_w$ method produces better fitness than its permanently constrained sibling (LM$_c$) in most cases, as the weaker approach lets model parameters vary freely unless mathematical coherence is being compromised. Both approaches generate the same outcome when the coherence of the bounds does break down and LM$_w$ effectively reverts to LM$_c$. Similarly, the PM approach may produce lower fit when it requires the Box-Cox transformation~\cite{box1964analysis,souza2017parametrized} of the regressand, to restore coherence of the bounds. We provide such an example in Appendix (see Fig.~\ref{fig:box_cox_PM}) showing relatively low fit of the PM approach in this case. Generally, both the CM and CIM methods yield poor model-fit.

In addition to generally comparing model fitness between methods, we also explored the impact of the key features (i.e., center and range) of the regressor on those of the regressand, considering an extreme case where the regressor is quasi certain against an uncertain regressand. In this instance, only the MinMax and PM methods produce desirable results, though we highlight that MinMax may break the coherence.

On the whole, the PM and LM models (particularly its weakly-constrained form, LM$_w$) appear to be the most widely suitable estimation approaches. Figure~\ref{fig:flowchart} provides a flowchart to summarise the findings presented. The figure provides a recommendation on suitability of regression approaches in respect to different features of the intervals in given data sets.

Lastly, as for computational efficiency, the LM$_c$ generally requires more execution time than LM$_w$ as it generically restricts parameters. Similarly, the PM approach involves extra execution time when it applies the Box-Cox transformation to restore consistency.

\begin{figure}
\centering
\includegraphics[width=1\columnwidth]{figure10.png}
\caption{A flowchart presenting recommendations on regression methods in respect to different interval-valued data sets.}
\label{fig:flowchart}
\end{figure}

\section{Conclusions}
\label{sec:conclusion}
This paper has carried out a comprehensive review and analysis of existing linear regression models, exploring their behaviour in respect to key features of both (systematically designed) synthetic and real-world interval-valued data sets, emphasising key features, such as whether they are puffy, skinny, scattered, densely placed or overlapping intervals. It has also investigated how these models do or do not maintain mathematical coherence, and the impact of ensuring coherence on model-fitness. It has proposed two constrained variants of the LM method~\cite{sun2015linear} to enhance its general suitability for real-world applications by ensuring mathematical coherence of interval bounds.

Going beyond the review and analysis of methods, the paper has discussed  the potential of intervals in regression and AI, in particular highlighting the value of mapping not only \emph{position} (as in the numeric case), but also range, for example to infer how uncertainty about a regressor impacts uncertainty or position of a regressand. To support adoption and facilitate interpretation of interval-valued regression results, the paper has also introduced intuitive visualization approaches, by means of 3D surface plots.

Finally, the paper has provided initial recommendations on which regression models are best suited to which type of interval-valued data sets, based on their features as outlined above. This is complemented by generic insights, such as generally better model prediction observed for interval sets with narrow range or condensed interval placement.
Analyses have shown that the PM~\cite{souza2017parametrized} method and the LM~\cite{sun2015linear} method (with its weaker set of restrictions, i.e., LM$_w$) are the most suitable regression approaches in many cases -- including where intervals are skinny, puffy or a mix thereof. The constrained variant of the LM model~\cite{sun2015linear}, LM$_c$ also performs equally well in many cases. However, when intervals are all densely placed and overlapping, other existing methods, such as CRM~\cite{neto2008centre}, CCRM~\cite{neto2010constrained}, and MinMax~\cite{billard2002symbolic} may turn into suitable estimation approaches, offering both relatively simple implementation and computation. In such instances, it is desirable to verify if the CRM~\cite{neto2008centre} and MinMax~\cite{billard2002symbolic} methods maintain consistency of the bounds. Certainly, maintaining this coherence automatically, such as by the PM, LM, and CCRM methods, tends to require more execution time. The latter may be particularly relevant as the number of variables increases.

With increasing popularity of interval-valued data -- to capture uncertainty and true range in data -- being able to analyse and derive the important insights available from these data using AI becomes an increasingly important topic. The substantial body of work advancing regression for interval-valued data over the last ten years, and advances such as those put forward in this paper, provide the foundation for leveraging intervals more generally in AI. Substantial work remains, from further refining our capacity for communicating resulting insights to advancing our modelling techniques, such as to non-linear approaches.

Beyond advancing techniques, in the future, we aim to explore additional real-world cases to demonstrate and advance the understanding of the value of using intervals, rather than `crisp', numeric data. Specifically, we will revisit the nature of the information captured by intervals, whether disjunctive (such as in confidence intervals), or conjunctive (true ranges, such as the numbers between 2 and 6) -- and explore how different analysis techniques could, and whether they should, deal with such data differently.

\begin{table*
\centering
\caption{Synthetic datasets used for sensitivity analysis (Set-1 to Set-11)}
\begin{subtable}{.24\linewidth}
\centering
\caption{Set-1}
\resizebox{1\columnwidth}{!}{
\begin{Huge}
\begin{tabular}{c|cc|cc} \hline
Entry&$X^-$	&$X^+$&$Y^-$&$Y^+$ \\\hline
1& 0.6	&0.95 &1.0	&2.0\\
2& 1.4	&1.92 &1.01	&3.01\\
3& 2.4	&3.13 &1.02	&4.02\\
4& 3.6	&4.45 &1.03	&5.03\\
5& 5.0  &5.95 &1.04 &6.04\\\hline
&$\mu_{X^c}$=2.94 & $\mu_{X^w}$=0.68 & $\mu_{Y^c}$=2.52  & $\mu_{Y^w}$=3.0 \\
& $\sigma_{X^c}$=1.87 & $\sigma_{X^w}$=0.24 & $\sigma_{Y^c}$=0.81 &  $\sigma_{Y^w}$=1.58\\\hline
\end{tabular}
\end{Huge}}
\end{subtable}
\begin{subtable}{.24\linewidth}
\centering
\caption{Set-2}
\resizebox{1\columnwidth}{!}{
\begin{Huge}
\begin{tabular}{c|cc|cc} \hline
Entry&$X^-$	&$X^+$&$Y^-$&$Y^+$ \\\hline
1& 0.6	&0.65 &1.0	&2.0\\
2& 1.4	&1.44 &1.01	&3.01\\
3& 2.4	&2.46 &1.02	&4.02\\
4& 3.6	&3.61 &1.03	&5.03\\
5& 5.0  &5.03 &1.04 &6.04\\\hline
& $\mu_{X^c}$=2.62 & $\mu_{X^w}$=0.04 & $\mu_{Y^c}$=2.52 &$\mu_{Y^w}$=3.0 \\
& $\sigma_{X^c}$=1.74 & $\sigma_{X^w}$=0.02 & $\sigma_{Y^c}$=0.81 &  $\sigma_{Y^w}$=1.58\\\hline
\end{tabular}
\end{Huge}}
\end{subtable}
\begin{subtable}{.24\linewidth}
\centering
\caption{Set-3}
\resizebox{1\columnwidth}{!}{
\begin{Huge}
\begin{tabular}{c|cc|cc} \hline
Entry&$X^-$	&$X^+$&$Y^-$&$Y^+$ \\\hline
1& 4.4	&4.7 &3.8	&4.2\\
2& 3.4	&3.8 &5.8	&6.3\\
3& 4.8	&5.3 &7.5	&7.8\\
4& 6.4	&6.7 &5.4	&5.8\\
5& 8.1  &8.5 &6.8   &7.3\\\hline
& $\mu_{X^c}$=5.61 & $\mu_{X^w}$=0.38& $\mu_{Y^c}$=6.07 &$\mu_{Y^w}$=0.42 \\
& $\sigma_{X^c}$=1.84 &$\sigma_{X^w}$=0.08& $\sigma_{Y^c}$=1.41&  $\sigma_{Y^w}$=0.08\\\hline
\end{tabular}
\end{Huge}}
\end{subtable}
\begin{subtable}{.24\linewidth}
\centering
\caption{Set-4}
\resizebox{1\columnwidth}{!}{
\begin{Huge}
\begin{tabular}{c|ll|ll} \hline
Entry&$X^-$	&$X^+$&$Y^-$&$Y^+$ \\\hline
1& 3.7&	4.0 & 5.5&	5.9\\
2& 3.4&	3.8 & 5.8&	6.3\\
3& 4.8&	5.3 & 7.5&	7.8\\
4& 7.9&	8.2 & 6.5&	6.9\\
5& 8.1&	8.5 & 6.8&	7.3\\\hline
& $\mu_{X^c}$=5.77 & $\mu_{X^w}$=0.38& $\mu_{Y^c}$=6.63 &$\mu_{Y^w}$=0.42 \\
& $\sigma_{X^c}$=2.26 &$\sigma_{X^w}$=0.08& $\sigma_{Y^c}$=0.78 &  $\sigma_{Y^w}$=0.08\\\hline
\end{tabular}
\end{Huge}}
\end{subtable}\vspace{1mm}\\%
\begin{subtable}{.24\linewidth}
\centering
\caption{Set-5}
\resizebox{1\columnwidth}{!}{
\begin{Huge}
\begin{tabular}{c|ll|ll} \hline
Entry&$X^-$	&$X^+$&$Y^-$&$Y^+$ \\\hline
1& 5.4&	5.7&	6.4&	6.8\\
2& 5.6&	6.0&	6.1&	6.6\\
3& 5.3&	5.8&	6.6&	6.9\\
4& 5.5&	5.8&	6.6&	7.0\\
5& 5.7&	6.1&	6.3&	6.8\\\hline
& $\mu_{X^c}$=5.69 & $\mu_{X^w}$=0.38& $\mu_{Y^c}$=6.61 &$\mu_{Y^w}$=0.42 \\
& $\sigma_{X^c}$=0.16 &$\sigma_{X^w}$=0.08& $\sigma_{Y^c}$=0.18 &  $\sigma_{Y^w}$=0.08\\\hline
\end{tabular}
\end{Huge}}
\end{subtable}
\begin{subtable}{.24\linewidth}
\centering
\caption{Set-6}
\resizebox{1\columnwidth}{!}{
\begin{Huge}
\begin{tabular}{c|ll|ll} \hline
Entry&$X^-$	&$X^+$&$Y^-$&$Y^+$ \\\hline
1& 1.2	&2.4	&4.8	&6.8\\
2& 2.6	&4.2	&3.6	&5.8\\
3& 3.8	&5.6	&7.6	&9.5\\
4& 5.0	&7.8	&2.6	&4.5\\
5& 7.6	&9.4	&5.0	&7.6\\\hline
& $\mu_{X^c}$=4.96 & $\mu_{X^w}$=1.84& $\mu_{Y^c}$=5.78 &$\mu_{Y^w}$=2.12 \\
& $\sigma_{X^c}$=2.60 &$\sigma_{X^w}$=0.59& $\sigma_{Y^c}$=1.88 &  $\sigma_{Y^w}$=0.29\\\hline
\end{tabular}
\end{Huge}}
\end{subtable}
\begin{subtable}{.24\linewidth}
\centering
\caption{Set-7}
\resizebox{1\columnwidth}{!}{
\begin{Huge}
\begin{tabular}{c|ll|ll} \hline
Entry&$X^-$	&$X^+$&$Y^-$&$Y^+$ \\\hline
1& 1.2&	2.4&	4.8&	6.8\\
2& 2.0&	3.6&	3.3&	5.5\\
3& 3.8&	5.6&	7.6&	9.5\\
4& 5.8&	8.6&	4.2&	6.1\\
5& 7.4&	9.2&	4.8&	7.4\\\hline
& $\mu_{X^c}$=4.96 & $\mu_{X^w}$=1.84& $\mu_{Y^c}$=6.00 &$\mu_{Y^w}$=2.12 \\
& $\sigma_{X^c}$=2.78 &$\sigma_{X^w}$=0.59& $\sigma_{Y^c}$=1.57 &  $\sigma_{Y^w}$=0.29\\\hline
\end{tabular}
\end{Huge}}
\end{subtable}%
\begin{subtable}{.24\linewidth}
\centering
\caption{Set-8}
\resizebox{1\columnwidth}{!}{
\begin{Huge}
\begin{tabular}{c|cc|cc} \hline
Entry&$X^-$	&$X^+$&$Y^-$&$Y^+$ \\\hline
1& 5.3&	6.5&	6.2&	8.2\\
2& 5.5&	7.1&	5.2&	7.4\\
3& 4.9&	6.7&	6.5&	8.4\\
4& 5.1&	7.9&	5.7&	7.6\\
5& 6.2&	8.0&	5.9&	8.5\\\hline
& $\mu_{X^c}$=6.32 & $\mu_{X^w}$=1.84& $\mu_{Y^c}$=6.96 &$\mu_{Y^w}$=2.12 \\
& $\sigma_{X^c}$=0.52 &$\sigma_{X^w}$=0.59& $\sigma_{Y^c}$=0.47 &  $\sigma_{Y^w}$=0.29\\\hline
\end{tabular}
\end{Huge}}
\end{subtable}\vspace{1mm}\\
\begin{subtable}{.24\linewidth}
\centering
\caption{Set-9}
\resizebox{1\columnwidth}{!}{
\begin{Huge}
\begin{tabular}{c|ll|ll} \hline
Entry&$X^-$	&$X^+$&$Y^-$&$Y^+$ \\\hline
1& 2.8&	5.4&	6.4&	9.2\\
2& 3.9&	5.1&	4.4&	5.8\\
3& 1.1&	3.2&	2.4&	4.6\\
4& 5.3&	5.6&	1.8&	2.3\\
5& 7.7&	9.3&	3.2&	5.0\\\hline
& $\mu_{X^c}$=4.94 & $\mu_{X^w}$=1.56 & $\mu_{Y^c}$=4.51 & $\mu_{Y^w}$=1.74 \\
& $\sigma_{X^c}$=2.32 & $\sigma_{X^w}$=0.88 & $\sigma_{Y^c}$=2.15&$\sigma_{Y^w}$=0.86\\\hline
\end{tabular}
\end{Huge}}
\end{subtable
\begin{subtable}{.24\linewidth}
\centering
\caption{Set-10}
\resizebox{1\columnwidth}{!}{
\begin{Huge}
\begin{tabular}{c|ll|ll} \hline
Entry&$X^-$	&$X^+$&$Y^-$&$Y^+$ \\\hline
1& 4.8&	6.4&	5.4&	7.2\\
2& 4.8&	5.1&	4.0&	4.5\\
3& 4.2&	5.4&	4.4&	5.8\\
4& 5.0&	7.6&	3.6&	6.4\\
5& 6.1&	8.2&	4.8&	7.0\\\hline
& $\mu_{X^c}$=5.76 & $\mu_{X^w}$=1.56& $\mu_{Y^c}$=5.31 &$\mu_{Y^w}$=1.74 \\
& $\sigma_{X^c}$=0.98 &$\sigma_{X^w}$=0.88& $\sigma_{Y^c}$=0.80 &  $\sigma_{Y^w}$=0.86\\\hline
\end{tabular}
\end{Huge}}
\end{subtable}%
\begin{subtable}{.24\linewidth}
\centering
\caption{Set-11}
\resizebox{1\columnwidth}{!}{
\begin{Huge}
\begin{tabular}{c|ll|ll} \hline
Entry&$X^-$	&$X^+$&$Y^-$&$Y^+$ \\\hline
1& 6.1&	6.4&	4.7&	5.2\\
2& 5.6&	6.8&	4.3&	5.7\\
3& 5.4&	7.0&	4.1&	5.9\\
4& 5.2&	7.3&	3.9&	6.1\\
5& 5.0&	7.6&	3.6&	6.4\\\hline
& $\mu_{X^c}$=6.24 & $\mu_{X^w}$=1.56& $\mu_{Y^c}$=4.99 &$\mu_{Y^w}$=1.74 \\
& $\sigma_{X^c}$=0.04 & $\sigma_{X^w}$=0.88 & $\sigma_{Y^c}$=0.02 &  $\sigma_{Y^w}$=0.86\\\hline
\end{tabular}
\end{Huge}}
\end{subtable}%
\label{tab:dataset2}
\end{table*}

\section*{Acknowledgment}
This work was supported by the UK EPSRC’s Leveraging the Multi-Stakeholder Nature of Cyber Security grant, EP/P011918/1.

\bibliographystyle{IEEEtran}
\balance