# A class of bootstrap based residuals for compositional data

Gustavo H. A. Pereira1,2 Jianwen Cai2

###### Abstract

Regression models for compositional data are common in several areas of
knowledge. As in other classes of regression models, it is desirable to
perform diagnostic analysis in these models using residuals that are
approximately standard normally distributed. However, for regression models
for compositional data, there has not been any multivariate residual that
meets this requirement. In this work, we introduce a class of asymptotically
standard normally distributed residuals for compositional data based on
bootstrap. Monte Carlo simulation studies indicate that the distributions of
the residuals of this class are well approximated by the standard normal
distribution in small samples. An application to simulated data also suggests
that one of the residuals of the proposed class is better to identify model
misspecification than its competitors. Finally, the usefulness of the best
residual of the proposed class is illustrated through an application on sleep
stages. The class of residuals proposed here can also be used in other classes
of multivariate regression models.

11footnotetext: Department of Statistics, Federal University of São Carlos,
Rod. Washington Luís, km 235 - SP-310 - São Carlos, CEP 13565-905, Brazil.
Email<EMAIL_ADDRESS>Department of Biostatistics, Gillings
School of Global Public Health, University of North Carolina at Chapel Hill,
Chapel Hill, NC, USA.

Keywords: Bootstrap, compositional data, diagnostic analysis, quantile
residual.

##  1 Introduction

Compositional data are proportions of disjoint categories adding to one
(Hijazi and Jernigan, 2009). They are common in several areas such as medicine
(Gueorguieva et al., 2008), marketing (Morais et al., 2018), ecology
(Tsilimigras and Fodor, 2016), genetics (Fernandes et al., 2014), and
geosciences (Filzmoser et al., 2012). This kind of data can be modeled in
different ways, but it is common to model them using the Dirichlet regression
(Hijazi and Jernigan, 2009) that can be viewed as a generalization of beta
regression (Ferrari and Cribari-Neto, 2004). An overview of advances in the
analysis of compositional data can be found in Alenazi (2023).

When a regression model is fitted, it is important to check model adequacy and
to identify outliers and influential observations. To reach these goals, it is
of interest to use residuals whose distribution is well approximated by the
standard normal distribution. When the distribution of a residual is not well
approximated by the standard normal distribution, it is hard to interpret
residuals plots. In these cases, it is not unusual that they are also
nonidentically distributed and well-fitted models may be discarded (see
Application 3 of Pereira (2019)).

The quantile residual (Dunn and Smyth, 1996) is simple and can be used in many
regression models. It is asymptotically standard normally distributed if the
parameters are consistently estimated. In addition, several papers showed that
its distribution in small samples is well approximated by the standard normal
distribution for different regression models (Pereira, 2019; Lemonte and
Moreno-Arenas, 2019; Feng et al., 2020). Other works proposed modifications in
the quantile residual to make it more suitable for some kinds of response
variables (Scudilio and Pereira, 2020; Pereira et al., 2020; Andrade et al.,
2023).

For compositional data, the quantile residual can be calculated in two
different ways. One can obtain the univariate or marginal quantile residual
calculating it for each component of each observation (Hijazi, 2006). The
problem with this approach is that it is hard to use a univariate residual for
identifying outliers and for checking model adequacy in the analysis of
compositional data, because each observation has several residuals. The other
option is to calculate a single multivariate quantile residual for each
observation using the distribution function assumed for the compositional
response variable. However, the multivariate quantile residual may not be a
good measure of the discrepancy between the observed and the fitted value of
the model.

Several other residuals were proposed for compositional data. Gueorguieva et
al. (2008) describes three univariate residuals for compositional data, but
they have the same drawback of the marginal quantile residual. Gueorguieva et
al. (2008) also proposed two multivariate residuals. They are known as
composite residuals, because they are a function of univariate residuals.
However, these residuals assume only positive values and so their distribution
cannot be well approximated by the standard normal distribution.

In this work, we propose a class of multivariate bootstrap based residuals for
compositional data and prove that they are asymptotically standard normally
distributed. Simulation studies also suggest that their distribution are well
approximated by the standard normal distribution in small samples. This class
of residuals use bootstrap to estimate the distribution function of a
multivariate residual and apply the quantile function of the standard normal
distribution for obtaining an asymptotically standard normally distributed
residual.

The remainder of this paper is organized as follows. Section 2 presents the
regression model considered in this work and introduces a class of
multivariate bootstrap based residuals for compositional data. In Section 3,
Monte Carlo simulation studies are performed to study the properties of four
residuals in the proposed class. An application to real data and to simulated
data are presented in Section 4. Concluding remarks are provided in Section 5.

##  2 A class of bootstrap based residuals

The class of residuals proposed in this section can be used in any parametric
regression model for compositional response variable and independent
observations. In the simulation studies and in the applications of our work,
we consider the Dirichlet regression model, because it is the most used
parametric model for compositional data. Some recent contributions related to
this model include Morais et al. (2018), Ross (2020), Hanretty (2021), and
Melo et al. (2022). The Dirichlet regression model has been widely used in
applications of several areas. See some recent applications of the model, for
example, in Vieira and Salmon (2019), Zhao et al. (2019), Chen et al. (2020),
Breuninger et al. (2021), and Reichert et al. (2022).

###  2.1 Dirichlet regression model

There are two parameterizations of the Dirichlet distribution and they are
known as common and alternative parameterization (Morais et al., 2018). In the
former we have $k$ shape parameters, where $k$ is the number of components of
the compositional random variable. Here we consider the alternative
parametrization in which the parameters are more interpretable, and so it is
more suitable for fitting regression models.

Let $w$ be a continuous compositional random variable, i.e.,
$w=(w_{1},w_{2},\ldots,w_{k})^{\top}$, $0<w_{j}<1$, $\sum_{j=1}^{k}w_{j}=1$.
The compositional random variable $w$ is Dirichlet distributed with parameters
$(\mu_{1},\mu_{2},\ldots,\mu_{k},\phi)^{\top},\mu_{1}=1-\sum_{j=2}^{k}\mu_{j}$,
$0<\mu_{j}<1,\forall j\in\\{1,2,\ldots,k\\}$, and $\phi>0$, if its probability
density function is given by (Da-Silva and Rodrigues, 2015)

$f(w_{1},\ldots,w_{k};\mu_{1},\ldots,\mu_{k},\phi)=\frac{\Gamma(\phi)}{\prod_{j=1}^{k}\Gamma(\phi\mu_{j})}\prod_{j=1}^{k}w_{j}^{\phi\mu_{j}-1}.$
(1)

In this parameterization, $\text{E}(w_{j})=\mu_{j}$ and
$\text{Var}(w_{j})=\frac{\mu_{j}(1-\mu_{j})}{1+\phi}$. Therefore the parameter
$\mu_{j}$ is the mean of the component $j$ of the compositional random
variable and $\phi$ is a precision parameter.

Let $y_{1},y_{2},\ldots,y_{n},y_{i}=(y_{i1},y_{i2},\ldots,y_{ik})^{\top}$ be
Dirichlet distributed random variables as defined in (1). The systematic
components of the Dirichlet logistic regression model are given by

$\begin{array}[]{lll}\log\left(\frac{\mu_{ij}}{\mu_{i1}}\right)&=&x_{ij1}\beta_{j1}+x_{ij2}\beta_{j2}+\ldots+x_{ijp_{j}}\beta_{jp_{j}},\,2\leq
j\leq k,\\\
\log(\phi_{i})&=&d_{i1}\gamma_{1}+d_{i2}\gamma_{2}+\ldots+d_{ip_{\phi}}\gamma_{p_{\phi}},\end{array}$
(2)

where
($x_{ij1},x_{ij2},\ldots,x_{ijp_{j}},d_{i1},d_{i2},\ldots,d_{ip_{\phi}})^{\top}$
are constants that represent the values of the explanatory variables and
$(\beta_{j1},\beta_{j2},\ldots,\beta_{jp_{j}},\gamma_{1},\gamma_{2},\ldots,\gamma_{p_{\phi}})^{\top}$
are unknown parameters. Note that model (2) uses the logit link function as is
used in the multinomial logistic regression (Hosmer Jr et al., 2013). As in
regression models with multinomial response (Das and Mukhopadhyay, 2014; Li et
al., 2022), other link functions can be used in regression models with
Dirichlet response. However, it is convenient to use the logit link function
for the equations related to the mean vector of the response variable and the
logarithmic link function for the dispersion parameter, because they lead to
interpretable parameters.

The parameters of model (2) can be estimated by maximum likelihood and the
model can be fitted using the package DirichletReg (Maier, 2014) of software
R. We used this package in the simulation studies and in the applications
presented in Sections 3 and 4, respectively.

###  2.2 A class of residuals for compositional data

For regression models for univariate responses and independent observations,
the quantile residual applies the inverse of the distribution function of the
standard normal distribution on the distribution function of the response
variable to obtain an asymptotically standard normally distributed residual.
If we do the same for compositional data, using the multivariate distribution
function of the response variable, we obtain a residual that is not
necessarily a good measure of the discrepancy between an observed and a fitted
value of the model. This fact is a consequence of the multivariate nature of
the compositional response variable.

To obtain an asymptotically standard normally distributed residual that is
also a good measure of the discrepancy between an observed and a fitted value
of the model, we modify the quantile residual. We change the distribution
function of the response variable by the distribution function of a
multivariate residual for compositional data. As the distribution of the
multivariate residuals for compositional data are unknown, we use bootstrap
(Efron and Tibshirani, 1994) to estimate the mentioned distribution function.
Therefore, we propose a class of residuals that can be obtained using the
following steps.

  1. 1.

For each observation, calculate $l_{i}=h_{i}r_{i}$, where $h_{i}$ is 1 or -1
(a sign function) and $r_{i}$ is a multivariate residual for compositional
data. In the end of this section, we suggest two options for $h_{i}$ and
$r_{i}$.

  2. 2.

Generate $B$ bootstrap samples of $y=(y_{1},y_{2},\ldots,y_{n})^{\top}$ using
the fitted model and denote each of them as
$y^{\\{b\\}}=(y_{1}^{\\{b\\}},y_{2}^{\\{b\\}},\ldots,y_{n}^{\\{b\\}})^{\top}$.

  3. 3.

For each bootstrap sample, fit the compositional regression model used for the
original data considering the same covariates of the original model.

  4. 4.

For each bootstrap sample, obtain
$l_{i}^{\\{b\\}}=h_{i}^{\\{b\\}}r_{i}^{\\{b\\}}$, where $h_{i}^{\\{b\\}}$ and
$r_{i}^{\\{b\\}}$ are $h_{i}$ and $r_{i}$, respectively, for the
$b{\text{th}}$ bootstrap sample.

  5. 5.

For each observation, obtain $a_{i}=\\#\\{l_{i}^{\\{b\\}}<l_{i}\\}$.

  6. 6.

For each observation, randomly generate $u_{i}$ from a uniform distribution in
the interval $(a_{i}/(B+1),(a_{i}+1)/(B+1))$.

  7. 7.

A residual of the proposed class is given by $s_{i}=\Phi^{-1}(u_{i})$, where
$\Phi$ is the cumulative distribution function of the standard normal
distribution.

Several observations related to the described procedure are necessary. The
first ones refer to step 1. As we can use as $r_{i}$ any multivariate residual
for compositional data and we can use different sign functions $h_{i}$,
$s_{i}$ is a class of residuals, because for each different kind of $r_{i}$
and $h_{i}$, we have a different kind of $s_{i}$.

In compositional regression models, to define a multivariate composite
residual, usually a univariate residual is chosen and the multivariate
residual is defined as the sum of the squared univariate residuals, across the
$k$ components of each observation (Gueorguieva et al., 2008). We use this way
of obtaining a multivariate residual, and also an alternative way in which we
replace the quadratic function by the absolute value function. For the
univariate residual, we consider the quantile residual, also known as pseudo
residual in compositional data (Hijazi, 2006), because it is suggested as the
best option to perform diagnostic analysis when the support of the response
variable is the unit interval (Pereira, 2019; Lemonte and Moreno-Arenas,
2019). The univariate quantile residual for the Dirichlet regression model is
defined as $q_{ij}=\Phi^{-1}[F(y_{ij};\hat{\mu}_{ij},\hat{\phi_{i}})]$, where
$\Phi(.)$ and $F(.)$ are the distribution function of the standard normal
distribution and of the beta distribution, respectively, and $\hat{\mu}_{ij}$
and $\hat{\phi}_{i}$ are the maximum likelihood estimator of $\mu_{ij}$ and
$\phi_{i}$, respectively. We denote the multivariate residual $r_{i}$ as
$r_{i}^{\\{a\\}}$ and $r_{i}^{\\{q\\}}$, when the absolute value function and
the quadratic function are used to calculate $r_{i}$, respectively.

Another observation is related to the sign function. It is important to use a
sign function $h_{i}$, because the possible choices for $r_{i}$ assume only
positive values and the best fitted observations have the lowest values of
$r_{i}$. Therefore, if the sign function was not used, $s_{i}$ would have a
high absolute value for the best fitted observations, because of the usage of
the function $\Phi(.)$ in the step 7 of the procedure to obtain the residual.
This is not an adequate behavior for a residual.

We consider in this work two different definitions for the sign function
$h_{i}$. For both definitions, we use an auxiliary random variable $m_{i}$ and
$h_{i}=\text{sign}\\{q_{im_{i}}\\}$. In the first, for each observation, we
consider the sign of the univariate quantile residual from the component that
has the largest distance between the observed and the fitted value (largest
absolute value of the univariate quantile residual). Therefore, in the first
definition of $h_{i}$ denoted by $h_{i}^{\\{1\\}}$,
$m_{i}^{\\{1\\}}=\underset{j\in\\{1,2,\ldots,k\\}}{\operatorname*{arg\,max}}|q_{ij}|$.
In the other definition, for all observations, we consider the sign of the
univariate quantile residual from the component that has the largest average
value of the response variable. Consequently, in this definition,
$m_{i}^{\\{2\\}}=\underset{j\in\\{1,2,\ldots,k\\}}{\operatorname*{arg\,max}}\sum_{i=1}^{n}y_{ij}/n$.

The last observation is related to steps 6 and 7. It is important to consider
$s_{i}=\Phi^{-1}(u_{i})$ and not obtaining $s_{i}$ as a direct function of
$a_{i}$ for two reasons. If we use $s_{i}$ as a direct function of $a_{i}$, we
could not obtain $s_{i}$ when $a_{i}$ equals to 0 or 1. The second reason is
to make $s_{i}$ continuous and allow it to converge in distribution to the
standard normal distribution using a finite $B$. Note that, as we usually use
a large $B$, the value of $s_{i}$ is not considerably affected by the
randomization made in step 6. This strategy is similar to that used to obtain
the randomized quantile residual for discrete data (Dunn and Smyth, 1996).

Considering that we suggested two possible definitions for $r_{i}$ and for
$h_{i}$, we consider in this work four different residuals in the proposed
class: $s_{i}^{\\{a1\\}}$, $s_{i}^{\\{q1\\}}$, $s_{i}^{\\{a2\\}}$, and
$s_{i}^{\\{q2\\}}$. Table 1 summarizes the features of this four residuals.

Table 1: Features of the residuals of the proposed class. Residual | Multivariate function used | Sign function
---|---|---
$s_{i}^{\\{a1\\}}$ | absolute value |  $h_{i}=\text{sign}\\{q_{im_{i}^{\\{1\\}}}\\}$, $m_{i}^{\\{1\\}}=\underset{j\in\\{1,2,\ldots,k\\}}{\operatorname*{arg\,max}}|q_{ij}|$
$s_{i}^{\\{q1\\}}$ | quadratic |  $h_{i}=\text{sign}\\{q_{im_{i}^{\\{1\\}}}\\}$, $m_{i}^{\\{1\\}}=\underset{j\in\\{1,2,\ldots,k\\}}{\operatorname*{arg\,max}}|q_{ij}|$
$s_{i}^{\\{a2\\}}$ | absolute value |  $h_{i}=\text{sign}\\{q_{im_{i}^{\\{2\\}}}\\}$, $m_{i}^{\\{2\\}}=\underset{j\in\\{1,2,\ldots,k\\}}{\operatorname*{arg\,max}}\sum_{i=1}^{n}y_{ij}/n$
$s_{i}^{\\{q2\\}}$ | quadratic |  $h_{i}=\text{sign}\\{q_{im_{i}^{\\{2\\}}}\\}$, $m_{i}^{\\{2\\}}=\underset{j\in\\{1,2,\ldots,k\\}}{\operatorname*{arg\,max}}\sum_{i=1}^{n}y_{ij}/n$

As mentioned before, the following theorem holds for the residuals of the
proposed class. It is proved in Appendix.

###### Theorem 2.2.1

Let $y_{1},...,y_{n}$ be compositional random variables, where each $y_{i}$
has vector of parameters $\tau_{i}$. If $\tau_{i}$ is consistently estimated,
under the true model,

$s_{i}\xrightarrow{D}N(0;1).$

Theorem 2.2.1 guarantees that the proposed residuals are approximately
standard normally distributed for large sample sizes. However, it is important
to study if this approximation can be made in small samples and to evaluate
their capacity to detect model misspecification. We perform these tasks in
Sections 3 and 4.1, respectively.

###  2.3 Other residuals for compositional data

The best known multivariate residual for compositional data is the composite
residual proposed by Gueorguieva et al. (2008). It was used for example by Yoo
et al. (2022) and can be calculated in the package DirichletReg of the
software R. This multivariate residual is a sum of the squared univariate
Pearson residuals across the $k$ components of each observation. In the
parameterization of the Dirichlet distribution used in this work, the
composite residual is given by

$r_{i}^{\\{q1\\}}=\sum_{j=1}^{k}\left(\frac{y_{ij}-\hat{\mu}_{ij}}{\sqrt{\hat{\phi}_{i}}}\right)^{2}.$
(3)

We also consider a residual similar to (3), but using the quantile residual as
the univariate residual. The composite quantile residual is given by

$r_{i}^{\\{q2\\}}=\sum_{j=1}^{k}q_{ij}^{2}.$ (4)

Residual (3) and (4) assume only positive values and so they are not
approximately standard normally distributed.

##  3 Simulation studies

Monte Carlo simulation studies were performed to study if the distribution of
the proposed residuals can be well approximated by the standard normal
distribution in small samples. For all of them, we considered two values for
the sample size (20 and 50), $k=3$, $B=1000$ and 2000 Monte Carlo replicates.
In each of the scenarios evaluated, we generated two covariates,
$x_{i22}=x_{i32}=d_{i2}$, $x_{i23}=x_{i33}=d_{i3}$ and
$x_{i21}=x_{i31}=d_{i1}=1$ for all $i$. The predictor variables were kept
fixed across the 2000 replicates. All simulations were performed using the
software R and the package DirichletReg. Given the goal of these simulation
studies, we did not include residuals (3) and (4) in these analysis.

In scenario 1a, the covariates were generated from the standard uniform
distribution and the value of the parameters were chosen in a way that the
mean of $\mu_{i1}$, $\mu_{i2}$, and $\mu_{i3}$ across the $n$ observations are
the same and the variance of the response variable is high ($\phi$ close to
20). In this scenario,
$\beta_{21}=-0.3,\beta_{22}=1.0,\beta_{23}=-0.5,\beta_{31}=-0.3,\beta_{32}=-0.5,\beta_{33}=1.0,\gamma_{1}=3.0,$
and $\gamma_{2}=\gamma_{3}=0$.

In scenarios 2a to 5a, we changed a single feature of the first one. In
scenarios 2a and 3a, we changed the value of some parameters to make the mean
of the components of the response variable far and very far from each other,
respectively. In scenario 4a, we considered that $\phi$ is a function of the
covariates using $\gamma_{2}=0.5$ and $\gamma_{3}=-0.5$. In the following
scenario, the first and second covariates were generated from the Bernoulli
distribution with parameter 0.5 and Gamma distribution with parameters 3 and 6
(mean and variance as the standard uniform distribution), respectively.
Scenarios 1b to 5b are similar to scenarios 1a to 5a, respectively, with
$\gamma_{1}=4.6$ ($\phi$ close to 100 that corresponds to a low variance of
the response variable). Table 2 summarizes the considered scenarios.

Table 2: Description of the scenarios. Scenarios |  | Mean of $(\mu_{1},\mu_{2},\mu_{3})$ |  | Covariates | $\phi$
---|---|---|---|---|---
1a, 1b |  | 0.333 | 0.333 | 0.333 |  | Uniform | Constant
2a, 2b |  | 0.290 | 0.151 | 0.559 |  | Uniform | Constant
3a, 3b |  | 0.308 | 0.049 | 0.643 |  | Uniform | Constant
4a, 4b |  | 0.333 | 0.333 | 0.333 |  | Uniform | Variable
5a, 5b |  | 0.333 | 0.333 | 0.333 |  | Bernoulli/Gamma | Constant

Table 3 presents the simulation results for scenario 1a and $n=20$. The table
presents the sample mean, variance, skewness, and kurtosis coefficients for
$s_{i}^{\\{a1\\}}$, $s_{i}^{\\{q1\\}}$, $s_{i}^{\\{a2\\}}$, and
$s_{i}^{\\{q2\\}}$ for each of the $n$ observations across the $2000$
replicates. The table also contains the value of the Anderson–Darling (AD)
statistic (Stephens, 1986) that can be used to test whether each residual is
standard normally distributed. The mean and the variance of residuals are very
close to zero and one, respectively, for all observations and proposed
residuals. The skewness coefficient is also very close to zero for almost all
observations and residuals and not so far from zero for the remaining
observations. The kurtosis coefficients of the residuals are farther for the
corresponding value of the standard normal distribution than the other
measures but there are no observations that have kurtosis coefficient very far
from 3 for any residual. As a consequence of these results, the values of the
AD statistics are all small and the distribution of the four proposed residual
seems to be well approximated by the standard normal distribution in this
scenario with $n=20$. The quality of this approximation seems to be similar
for the four residuals since the mean of the AD statistics are very close to
each other.

Table 4 presents the average values for the mean, variance, skewness, kurtosis
coefficient, and AD Statistics for all proposed residuals in the ten studied
scenarios with $n=20$. The distribution of the four proposed residuals seems
to be well approximated by the standard normal distribution in all scenarios
when $n=20$, since the average values of all measures for the four proposed
residuals are very close to theoretical values of this distribution. The
average values of the AD statistics are small for the four residuals in the
ten scenarios suggesting that this approximation is not affected by the value
of the parameters related to the mean of the response variable, nor by the
variance of the response and nor by the distribution of the covariates. The
quality of the standard normality approximation seems to be similar for the
four residuals in all scenarios. Results for $n=50$ are similar and not
included here for brevity.

##  4 Applications

We present two applications. The first application uses simulated data and its
goal is to compare the ability of the residuals to detect model
misspecification. We also analyzed data on sleep stages to illustrate the
usefulness of the best residual based on the simulation results in the first
application.

###  4.1 Application to simulated data

For evaluating if a parametric regression model properly fits a certain
response variable, a normal probability plot with simulated envelope
(Atkinson, 1981) is commonly used. If all residuals are inside the envelope or
if only few points are outside the envelope and they are all close to one of
the bounds of the envelope, there is no evidence that the model is inadequate.

Table 3: Simulation results Scenario 1a - $n=20$.

$i$ | Mean |  | Variance |  | Skewness |  | Kurtosis |  | A-D Statistic
---|---|---|---|---|---|---|---|---|---
| $s_{i}^{\\{a1\\}}$ | $s_{i}^{\\{q1\\}}$ | $s_{i}^{\\{a2\\}}$ | $s_{i}^{\\{q2\\}}$ |  | $s_{i}^{\\{a1\\}}$ | $s_{i}^{\\{q1\\}}$ | $s_{i}^{\\{a2\\}}$ | $s_{i}^{\\{q2\\}}$ |  | $s_{i}^{\\{a1\\}}$ | $s_{i}^{\\{q1\\}}$ | $s_{i}^{\\{a2\\}}$ | $s_{i}^{\\{q2\\}}$ |  | $s_{i}^{\\{a1\\}}$ | $s_{i}^{\\{q1\\}}$ | $s_{i}^{\\{a2\\}}$ | $s_{i}^{\\{q2\\}}$ |  | $s_{i}^{\\{a1\\}}$ | $s_{i}^{\\{q1\\}}$ | $s_{i}^{\\{a2\\}}$ | $s_{i}^{\\{q2\\}}$
1 | 0.01 | 0.01 | 0.01 | 0.01 |  | 0.94 | 0.94 | 0.95 | 0.94 |  | 0.03 | 0.06 | 0.00 | -0.03 |  | 2.96 | 3.09 | 3.02 | 3.01 |  | 0.88 | 0.94 | 0.76 | 0.82
2 | -0.01 | -0.02 | 0.02 | 0.02 |  | 1.06 | 1.06 | 1.06 | 1.06 |  | 0.13 | 0.10 | -0.01 | -0.01 |  | 2.96 | 2.88 | 2.98 | 2.93 |  | 1.43 | 1.45 | 1.23 | 1.21
3 | 0.01 | 0.01 | 0.04 | 0.04 |  | 1.05 | 1.05 | 1.06 | 1.06 |  | -0.07 | -0.07 | -0.02 | 0.01 |  | 3.09 | 3.16 | 3.20 | 3.22 |  | 0.72 | 0.83 | 1.63 | 1.77
4 | -0.02 | -0.02 | -0.01 | -0.01 |  | 0.97 | 0.98 | 0.97 | 0.98 |  | -0.06 | -0.04 | 0.03 | 0.03 |  | 2.91 | 2.90 | 2.90 | 2.90 |  | 0.52 | 0.49 | 0.46 | 0.40
5 | 0.02 | 0.02 | -0.01 | -0.01 |  | 0.98 | 0.96 | 0.98 | 0.97 |  | 0.11 | 0.10 | -0.04 | -0.03 |  | 3.09 | 2.99 | 3.13 | 3.05 |  | 1.19 | 1.16 | 0.74 | 0.75
6 | 0.00 | 0.00 | -0.01 | -0.01 |  | 0.99 | 0.99 | 0.99 | 0.99 |  | 0.05 | 0.05 | -0.03 | -0.04 |  | 3.31 | 3.37 | 3.40 | 3.40 |  | 1.09 | 1.21 | 1.14 | 1.25
7 | 0.01 | 0.01 | -0.01 | -0.01 |  | 1.04 | 1.04 | 1.04 | 1.04 |  | 0.04 | 0.00 | -0.17 | -0.14 |  | 3.13 | 3.17 | 3.10 | 3.09 |  | 0.67 | 0.75 | 0.98 | 0.92
8 | 0.00 | 0.00 | -0.03 | -0.04 |  | 1.02 | 1.02 | 1.02 | 1.02 |  | 0.03 | 0.03 | 0.01 | -0.03 |  | 3.02 | 2.99 | 3.03 | 2.98 |  | 0.42 | 0.39 | 1.13 | 1.30
9 | 0.03 | 0.04 | -0.01 | -0.01 |  | 1.02 | 1.03 | 1.03 | 1.03 |  | 0.09 | 0.11 | 0.01 | 0.00 |  | 3.06 | 3.11 | 3.14 | 3.08 |  | 1.17 | 1.20 | 0.47 | 0.40
10 | 0.00 | -0.01 | 0.02 | 0.02 |  | 1.02 | 1.02 | 1.02 | 1.02 |  | 0.08 | 0.06 | -0.06 | -0.06 |  | 2.99 | 3.02 | 3.04 | 3.07 |  | 0.61 | 0.58 | 0.92 | 0.92
11 | -0.02 | -0.03 | 0.01 | 0.01 |  | 0.95 | 0.95 | 0.95 | 0.95 |  | 0.05 | 0.02 | -0.07 | -0.06 |  | 3.00 | 3.03 | 3.06 | 3.08 |  | 1.37 | 1.38 | 0.86 | 0.89
12 | -0.02 | -0.02 | 0.02 | 0.02 |  | 0.99 | 0.99 | 0.99 | 0.99 |  | 0.00 | -0.01 | -0.04 | -0.09 |  | 2.98 | 2.95 | 2.90 | 2.94 |  | 0.46 | 0.50 | 0.86 | 0.91
13 | 0.01 | 0.00 | 0.02 | 0.02 |  | 1.02 | 1.03 | 1.03 | 1.03 |  | 0.01 | -0.01 | 0.01 | -0.01 |  | 2.94 | 2.96 | 2.92 | 2.94 |  | 0.71 | 0.63 | 0.88 | 0.86
14 | 0.04 | 0.03 | -0.03 | -0.03 |  | 0.99 | 0.99 | 1.00 | 1.01 |  | -0.03 | -0.03 | 0.02 | 0.01 |  | 3.01 | 3.03 | 3.05 | 3.11 |  | 1.85 | 1.90 | 1.04 | 1.01
15 | 0.03 | 0.03 | -0.04 | -0.04 |  | 1.00 | 1.00 | 1.01 | 1.00 |  | 0.08 | 0.07 | -0.01 | -0.02 |  | 2.91 | 2.89 | 3.02 | 2.95 |  | 0.95 | 0.91 | 1.39 | 1.44
16 | 0.01 | 0.01 | -0.01 | -0.01 |  | 0.97 | 0.98 | 0.98 | 0.98 |  | 0.07 | 0.03 | -0.03 | -0.01 |  | 2.82 | 2.84 | 2.89 | 2.86 |  | 0.30 | 0.21 | 0.30 | 0.23
17 | -0.02 | -0.03 | 0.02 | 0.01 |  | 0.99 | 0.98 | 0.98 | 0.98 |  | 0.13 | 0.09 | 0.03 | 0.00 |  | 3.26 | 3.01 | 3.02 | 3.03 |  | 1.13 | 1.08 | 0.81 | 0.79
18 | -0.02 | -0.02 | 0.00 | 0.00 |  | 1.03 | 1.04 | 1.03 | 1.03 |  | -0.05 | -0.02 | -0.04 | -0.03 |  | 2.95 | 2.93 | 2.87 | 2.94 |  | 0.70 | 0.65 | 0.41 | 0.43
19 | 0.02 | 0.02 | 0.00 | 0.00 |  | 0.94 | 0.95 | 0.96 | 0.96 |  | 0.01 | -0.03 | -0.05 | -0.07 |  | 2.94 | 3.04 | 3.01 | 3.07 |  | 1.01 | 0.95 | 0.54 | 0.53
20 | 0.00 | 0.00 | -0.04 | -0.04 |  | 1.02 | 1.02 | 1.02 | 1.02 |  | 0.07 | 0.05 | -0.03 | -0.02 |  | 2.90 | 2.89 | 2.97 | 2.94 |  | 0.67 | 0.62 | 1.76 | 1.95
Mean | 0.00 | 0.00 | 0.00 | 0.00 |  | 1.00 | 1.00 | 1.00 | 1.00 |  | 0.04 | 0.03 | -0.03 | -0.03 |  | 3.01 | 3.01 | 3.03 | 3.03 |  | 0.89 | 0.89 | 0.92 | 0.94
SD | 0.02 | 0.02 | 0.02 | 0.02 |  | 0.03 | 0.04 | 0.03 | 0.03 |  | 0.06 | 0.05 | 0.05 | 0.04 |  | 0.12 | 0.13 | 0.12 | 0.12 |  | 0.39 | 0.41 | 0.39 | 0.45

Table 4: Average distributional measures for the residuals - $n=20$. Scenarios |  | $\gamma_{1}=3.0$ |  | $\gamma_{1}=4.6$
---|---|---|---|---
|  | $s_{i}^{\\{a1\\}}$ | $s_{i}^{\\{q1\\}}$ | $s_{i}^{\\{a2\\}}$ | $s_{i}^{\\{q2\\}}$ |  | $s_{i}^{\\{a1\\}}$ | $s_{i}^{\\{q1\\}}$ | $s_{i}^{\\{a2\\}}$ | $s_{i}^{\\{q2\\}}$
| Mean | 0.00 | 0.00 | 0.00 | 0.00 |  | 0.00 | 0.00 | 0.00 | 0.00
| Variance | 1.00 | 1.00 | 1.00 | 1.00 |  | 1.00 | 1.00 | 1.00 | 1.00
1 | Skewness | 0.04 | 0.03 | -0.03 | -0.03 |  | 0.00 | -0.01 | 0.01 | 0.00
| Kurtosis | 3.01 | 3.01 | 3.03 | 3.03 |  | 2.98 | 3.01 | 2.96 | 2.99
| A-D Statistics | 0.89 | 0.89 | 0.92 | 0.94 |  | 0.90 | 0.87 | 0.92 | 0.92
| Mean | 0.00 | 0.00 | 0.00 | 0.00 |  | 0.00 | -0.01 | 0.00 | 0.00
| Variance | 1.00 | 1.00 | 1.00 | 1.00 |  | 1.00 | 1.00 | 1.00 | 1.00
2 | Skewness | 0.01 | 0.00 | 0.01 | 0.01 |  | 0.00 | 0.00 | -0.01 | -0.01
| Kurtosis | 2.99 | 3.00 | 3.00 | 2.99 |  | 2.99 | 2.99 | 3.01 | 3.01
| A-D Statistics | 0.91 | 0.90 | 0.98 | 0.99 |  | 0.99 | 0.99 | 0.80 | 0.80
| Mean | 0.00 | 0.00 | 0.00 | 0.00 |  | 0.01 | 0.01 | 0.00 | 0.00
| Variance | 1.00 | 1.00 | 1.00 | 1.00 |  | 1.00 | 1.00 | 1.00 | 1.00
3 | Skewness | -0.01 | -0.01 | 0.00 | 0.01 |  | 0.00 | 0.00 | 0.00 | 0.00
| Kurtosis | 2.99 | 3.00 | 2.97 | 2.97 |  | 3.01 | 3.02 | 2.99 | 3.01
| A-D Statistics | 1.24 | 1.27 | 0.81 | 0.79 |  | 0.84 | 0.83 | 0.86 | 0.81
| Mean | 0.01 | 0.01 | 0.00 | 0.00 |  | 0.01 | 0.01 | 0.00 | 0.00
| Variance | 0.99 | 1.00 | 1.00 | 1.00 |  | 1.00 | 1.00 | 1.00 | 1.00
4 | Skewness | 0.01 | 0.00 | -0.01 | 0.00 |  | 0.00 | 0.00 | 0.00 | 0.00
| Kurtosis | 2.95 | 2.98 | 2.99 | 3.00 |  | 2.99 | 3.01 | 2.97 | 2.99
| A-D Statistics | 1.46 | 1.46 | 1.22 | 1.18 |  | 0.91 | 0.89 | 1.17 | 1.11
| Mean | 0.01 | 0.01 | 0.00 | 0.00 |  | 0.00 | 0.00 | 0.00 | 0.00
| Variance | 1.00 | 1.00 | 1.00 | 1.00 |  | 1.00 | 1.00 | 1.00 | 1.00
5 | Skewness | 0.01 | 0.01 | -0.01 | -0.01 |  | 0.00 | 0.00 | -0.01 | -0.02
| Kurtosis | 3.00 | 3.01 | 3.03 | 3.01 |  | 2.97 | 2.97 | 2.97 | 2.97
| A-D Statistics | 0.65 | 0.66 | 0.90 | 0.91 |  | 1.08 | 1.05 | 0.88 | 0.85

There is no general rule in the literature to define how far a residual should
be to the closest bound of the envelope to lead to rejection of the model.
However, for a proper comparison of the ability of the residuals to detect
model misspecification, it is important to define a rule for model rejection.
Let $v$ be a constant for which, when the correct model is fitted, the
probability of at least one residual is outside the envelope and is more than
$v$ units away from the closest bound of the envelope is $5\%$. In practice,
with a given data set, $v$ will depend on the sample size $n$ and the
dispersion of the residuals.

Consider that $g$ datasets of size $n$ were generated from a certain
regression model. Consider also that for each of these datasets, the correct
model was fitted, a certain residual $r_{i}$ was calculated for the $n$
observations, and a normal probability plot with simulated envelope was made.
Let $e_{j}$ be zero if all residuals of the $j$th datasets are inside the
simulated envelope and be the distance from the residual that is farthest from
the envelope to the closest bound of the envelope. Using this definition, a
natural estimator of $v$ is given by $(e_{[0.95g]}+e_{[0.95g+1]})/2$, where
$e_{[j]}$ denotes the $j$th lowest value of $e_{j}$ across the $g$ datasets.

To estimate $v$ for the residuals considered in this work, we generated 500
datasets based on model (2) with $n=50$ and $k=3$. The two covariates were
generated from the binomial distribution with parameters $n$ and 0.5 and the
Dirichlet distributed response variable was generated according to the set 1
of Distributions presented on Table 5.

Table 5: Probability distributions used in the application with simulated data. $x_{i22}$ | $x_{i23}$ | Set 1 of distributions | Set 2 of distributions
---|---|---|---
1 | 0 | Dirichlet$(0.80,0.15,0.05,e^{4.6})$ | Dirichlet$(0.40,0.35,0.25,e^{4.6})$
0 | 1 | Dirichlet$(0.05,0.80,0.15,e^{4.6})$ | Dirichlet$(0.25,0.40,0.35,e^{4.6})$
1 | 1 | Dirichlet$(0.15,0.05,0.80,e^{4.6})$ | Dirichlet$(0.35,0.25,0.40,e^{4.6})$
0 | 0 | Dirichlet$(0.35,0.35,0.30,e^{4.6})$ | Dirichlet$(0.30,0.30,0.35,e^{4.6})$

For each of the datasets, we fitted the correct model, calculated the
residuals presented in Section 2, and made normal probability plots with
simulated envelope. As the residuals of the proposed class have a similar
variance as presented in Table 4, we used the same $v$ for these four
residuals. To estimate it, we used the plots for these four residuals and
denoted its natural estimator by $v^{\\{1\\}}$. Residuals (3) and (4) are a
sum of squared univariate residuals, and they have the same variances but are
different from those from the proposed class. Therefore, we used the same $v$
for these two residuals and denoted its estimator as $v^{\\{2\\}}$.

Using the 500 generated datasets, Figure 1 presents a histogram for the number
of points outside the simulated envelope and $v^{\\{j\\}}$ units or more away
from the closest envelope bounds, where $j=1$ for the residuals of the
proposed class and $j=2$ for $r_{i}^{\\{q1\\}}$ and $r_{i}^{\\{q2\\}}$. As
expected, the plots are similar and, for the six residuals, about $95\%$ of
datasets do not have any point outside the simulated envelope and
$v^{\\{j\\}}$ units or more away from the closest envelope bounds.

Figure 1: Histogram for the number of points outside the simulated envelope
and $v^{\\{j\\}}$ units or more away from the closest envelope bounds -
Correct model.

To compare the ability of the residuals to detect model misspecification, we
generated 500 additional datasets with $k=3$ and $n=50$. The covariates were
generated in the same way as in the previous datasets, and the response
variable was generated as a mixture of the set of Distributions 1 and 2 (Table
5) with probabilities 0.7 and 0.3, respectively. For each of these datasets,
we incorrectly fitted model (2), calculated the residuals presented in Section
2, and made normal probability plots with simulated envelope.

Figure 2 presents a histogram for the number of points outside the simulated
envelope and $v^{\\{j\\}}$ units or more away from the closest envelope bounds
when the incorrect model was fitted. The residual $s_{i}^{\\{a1\\}}$ of the
proposed class is the best to detect model misspecification since, for this
residual, only $11\%$ of the datasets do not have any point outside the
envelope and distant from the bounds. In addition, using $s_{i}^{\\{a1\\}}$,
more than $70\%$ of the datasets have at least 3 points outside the envelope
and far from the bounds. On the other hand, the other residuals could not
detect model misspecification for more than $20\%$ of the datasets.

Figure 2: Histogram for the number of points outside the simulated envelope
and $v^{\\{j\\}}$ units or more away from envelope bounds - Wrong model.

###  4.2 Application to data on sleep stages

The human sleep has four stages: N1, N2, N3 and REM. During the first two
stages, the person is in a light sleep (Ritmala-Castren et al., 2015). The
stage N3 corresponds a deep sleep and it is when muscles and tissues are
healed (Hussain et al., 2022). Finally, the stage REM is characterized by
rapid eye movements and fast breathing and this stage is responsible for
optimal emotional processing (Tempesta et al., 2018).

Considering the functions of the last two stages, it is important that people
spend a considerable percentage of their sleep on stages N3 and REM. Sleep is
considered as within normal parameters when the percentage time spent on the
stages N1, N2, N3, and REM are on the following ranges: $3-8\%$, $45-55\%$,
$15-20\%$, $20-25\%$, respectively (González-Naranjo et al., 2019).

The proportion of time spent on each sleep stage is a compositional variable.
However, papers usually study each of the stages separately considering four
independent proportions (Ohayon et al., 2017; Maski et al., 2021). For the
best of our knowledge, this paper is the first that considers the proportion
of time spent on each sleep stage as a compositional variable.

Veje et al. (2021) collected sleep data from 42 adults, in which 22 were
patients with Tick-borne encephalitis (TBE) and 20 controls. Their main goal
was to investigate if the TBE affects the quality of the sleep measured by a
questionnaire. Here, we use a Dirichlet regression model to study if the
proportion of time spent on each sleep stages changes with the total sleep
time (TST) and from patients with TBE to healthy adults. The sample means for
the proportion of time spent on N1, N2, N3 and REM stages are 0.134, 0.532,
0.152, and 0.182, respectively. These sample means are close for TBE adults
and control and the standard deviations are less than 0.1.

The Dirichlet regression model assumes that $P(y_{ij}=0)=0$ for all $i$ and
$j$. As a result, the model can not be fitted if there are one or more zeros
in the response variable. When there are many zeros, a possible solution is to
fit a zero adjusted Dirichlet regression model (Tsagris and Stewart, 2018).
However, the parameters of the discrete component of this model can not be
well estimated when there are few zeros. In this case, it is better to
transform data in a way that eliminates the zeros (Hijazi, 2011). When there
are very few zeros, it is reasonable to simply replacing the zeros by a small
value and then normalizing the imputed compositions (Aitchison, 1986).

Two patients of the data collected by Veje et al. (2021) do not spend time on
the N3 stage in their monitored night of sleep. We replaced their zero values
by 0.001 and normalized the proportions of these patients to sum 1.

We fitted a Dirichlet regression model for the transformed data containing the
mentioned covariates in all submodels. Using the log likelihood ratio test, we
concluded that there are no evidences that the parameters of the proportion of
time spent on each sleep stage changes from patients with TBE to healthy
adults (p = 0.7469). However, the mean vector (p = 0.0019) and the precision
parameter (p = 0.0030) of the response variable is a function of the total
sleep time.

As the residual $s_{i}^{\\{a1\\}}$ performed better in the analysis with
simulated data, we used it to conduct the diagnostic analysis. Figure 3
presents a normal probability plot with simulated envelope for
$s_{i}^{\\{a1\\}}$ considering the final model. The plot does not suggest
model misspecification of the Dirichlet regression model for the proportion of
time spent on each sleep stage with the total sleep time as a covariate of all
submodels.

Figure 3: Normal probability plot with simulated envelope for
$s_{i}^{\\{a1\\}}$ considering the final model for data on sleep stages

Table 6 presents the parameter estimates and standard errors for the final
model. Note that the estimates of the parameters associated with TST in the
second and third submodels are positive. Therefore, on average, an increase in
the total sleep time leads to a higher proportion of time spent on stages N3
and REM in relation to the N1 stage. These results suggest that adults have a
better quality of sleep if they spend enough time on sleeping.

Calculating the exponential of the parameter estimates (last column of Table
6), it is estimated that for every hour increase in the total sleep time, the
ratio between the mean of the proportion of time spent on the N3 stage of the
sleep and the mean of the proportion of time spent on the N1 stage increases
by 40.3%. Similarly, it is estimated that for every hour increase in the total
sleep time, the ratio between the mean of the proportion of time spent on the
REM stage of the sleep and the mean of the proportion of time spent on the N1
stage increases by 23.8%

Table 6: Parameter estimates and standard errors in the model for the proportion of time spent on each sleep stage. Submodel | Covariate | Estimate | Std. Error | Exp(estim)
---|---|---|---|---
$\mu_{N2}$ | Intercept | 0.5142 | 0.5949 | 1.672
| TST | 0.1314 | 0.0880 | 1.140
$\mu_{N3}$ | Intercept | -2.2371 | 0.7087 | 0.107
| TST | 0.3386 | 0.1030 | 1.403
$\mu_{REM}$ | Intercept | -1.0275 | 0.6688 | 0.358
| TST | 0.2136 | 0.0981 | 1.238
$\phi$ | Intercept | 0.8401 | 0.7433 | 2.317
| TST | 0.3400 | 0.1139 | 1.405

##  5 Concluding remarks

In this work, we proposed a class of asymptotically standard normally
distributed residuals for compositional data based on bootstrap. We
investigated the properties of four residuals of the proposed class in small
samples using Monte Carlo simulation studies and an application on simulated
data considering the Dirichlet regression model. We also presented an
application on the stages of the sleep that illustrates the usefulness of one
of the residuals of the proposed class.

To facilitate the diagnostic analysis of a regression model, it is desirable
to find a residual whose distribution is well approximated by the standard
normal distribution. To fill in the gap in the literature for compositional
data, we proposed a class of residuals and proved that the residuals of the
proposed class are asymptotically standard normally distributed, and hence
their distribution is close to the standard normal distribution when the
sample size is large. In addition, simulation studies suggested that the
distributions of the four investigated residual of the proposed class are very
close to the standard normal distribution even when $n=20$ for all considered
scenarios.

It is also important that a residual can detect model misspecification. An
application using simulated data suggested that the residual of the proposed
class denoted by $s_{i}^{\\{a1\\}}$ is good for identifying model
misspecification and better than the other residuals considered in this paper.

Another interesting feature of the proposed class of residuals is that their
expressions do not depend on the assumed distribution of the response
variable. The algebraic expression of the multivariate residuals based on the
deviance residual, for example, have to be obtained for each different
distribution assumed to the response variable.

Besides the usage for diagnostic analysis in regression models for
compositional data, the proposed class of residuals can also be used in other
classes of multivariate regression models. There is no approximately standard
normally distributed multivariate residual in the literature for multivariate
regression models. Therefore, our class of residuals can be a useful tool for
several classes of multivariate regression models. Future work can investigate
the properties of residuals of the proposed class in different multivariate
regression models, such as the multivariate normal regression model (Johnson
and Wichern, 2013), the multivariate elliptical regression model (Lemonte and
Patriota, 2011), and the multivariate generalized linear model (Schaid et al.,
2019).

##  6 Acknowledgment

This work was partially supported by São Paulo Research Foundation (FAPESP),
grant number 2020/16334-9.

## References

  * Aitchison (1986) Aitchison, J., 1986.  The Statistical Analysis of Compositional Data.  Springer.
  * Alenazi (2023) Alenazi, A., 2023.  A review of compositional data analysis and recent advances.  Communications in Statistics-Theory and Methods 52, 5535–5567.
  * Andrade et al. (2023) Andrade, A.C., Pereira, G.H., Artes, R., 2023.  The circular quantile residual.  Computational Statistics & Data Analysis 178, 107612.
  * Atkinson (1981) Atkinson, A.C., 1981.  Two graphical displays for outlying and influential observations in regression.  Biometrika 68, 13–20.
  * Breuninger et al. (2021) Breuninger, T.A., Wawro, N., Breuninger, J., Reitmeier, S., Clavel, T., Six-Merker, J., Pestoni, G., Rohrmann, S., Rathmann, W., Peters, A., et al., 2021.  Associations between habitual diet, metabolic disease, and the gut microbiota using latent dirichlet allocation.  Microbiome 9, 1–18.
  * Chen et al. (2020) Chen, B., Bernard, J.Y., Padmapriya, N., Ning, Y., Cai, S., Lança, C., Tan, K.H., Yap, F., Chong, Y.S., Shek, L., et al., 2020.  Associations between early-life screen viewing and 24 hour movement behaviours: findings from a longitudinal birth cohort study.  The Lancet Child & Adolescent Health 4, 201–209.
  * Da-Silva and Rodrigues (2015) Da-Silva, C.Q., Rodrigues, G.S., 2015.  Bayesian dynamic dirichlet models.  Communications in Statistics-Simulation and Computation 44, 787–818.
  * Das and Mukhopadhyay (2014) Das, I., Mukhopadhyay, S., 2014.  On generalized multinomial models and joint percentile estimation.  Journal of Statistical Planning and Inference 145, 190–203.
  * Dunn and Smyth (1996) Dunn, P.K., Smyth, G.K., 1996.  Randomized quantile residuals.  Journal of Computational and Graphical Statistics 5, 236–244.
  * Efron and Tibshirani (1994) Efron, B., Tibshirani, R.J., 1994.  An introduction to the bootstrap.  CRC press.
  * Feng et al. (2020) Feng, C., Li, L., Sadeghpour, A., 2020.  A comparison of residual diagnosis tools for diagnosing regression models for count data.  BMC Medical Research Methodology 20, 1–21.
  * Fernandes et al. (2014) Fernandes, A.D., Reid, J.N., Macklaim, J.M., McMurrough, T.A., Edgell, D.R., Gloor, G.B., 2014.  Unifying the analysis of high-throughput sequencing datasets: characterizing rna-seq, 16s rrna gene sequencing and selective growth experiments by compositional data analysis.  Microbiome 2, 1–13.
  * Ferrari and Cribari-Neto (2004) Ferrari, S., Cribari-Neto, F., 2004.  Beta regression for modelling rates and proportions.  Journal of applied statistics 31, 799–815.
  * Filzmoser et al. (2012) Filzmoser, P., Hron, K., Reimann, C., 2012.  Interpretation of multivariate outliers for compositional data.  Computers & Geosciences 39, 77–85.
  * González-Naranjo et al. (2019) González-Naranjo, J.E., Alfonso-Alfonso, M., Grass-Fernandez, D., Morales-Chacón, L.M., Pedroso-Ibáñez, I., Ricardo-De La Fe, Y., Padrón-Sánchez, A., 2019.  Analysis of sleep macrostructure in patients diagnosed with parkinson’s disease.  Behavioral Sciences 9, 6.
  * Gueorguieva et al. (2008) Gueorguieva, R., Rosenheck, R., Zelterman, D., 2008.  Dirichlet component regression and its applications to psychiatric data.  Computational statistics & data analysis 52, 5344–5355.
  * Hanretty (2021) Hanretty, C., 2021.  Forecasting multiparty by-elections using dirichlet regression.  International Journal of Forecasting 37, 1666–1676.
  * Hijazi (2011) Hijazi, R., 2011.  An em-algorithm based method to deal with rounded zeros in compositional data under dirichlet models, in: Proceedings of CoDaWork’11: 4th international workshop on Compositional Data Analysis, Egozcue, JJ, Tolosana-Delgado, R. and Ortego, MI (eds.) 2011, CIMNE. pp. 1–10.
  * Hijazi (2006) Hijazi, R.H., 2006.  Residuals and diagnostics in dirichlet regression.  ASA Proceedings of the General Methodology Section , 1190–1196.
  * Hijazi and Jernigan (2009) Hijazi, R.H., Jernigan, R.W., 2009.  Modelling compositional data using dirichlet regression models.  Journal of Applied Probability & Statistics 4, 77–91.
  * Hosmer Jr et al. (2013) Hosmer Jr, D.W., Lemeshow, S., Sturdivant, R.X., 2013.  Applied logistic regression.  3 ed., John Wiley & Sons.
  * Hussain et al. (2022) Hussain, I., Hossain, M.A., Jany, R., Bari, M.A., Uddin, M., Kamal, A.R.M., Ku, Y., Kim, J.S., 2022.  Quantitative evaluation of eeg-biomarkers for prediction of sleep stages.  Sensors 22, 3079.
  * Johnson and Wichern (2013) Johnson, R.A., Wichern, D.W., 2013.  Applied multivariate statistical analysis.  6 ed., Prentice hall Upper Saddle River, NJ.
  * Lemonte and Moreno-Arenas (2019) Lemonte, A.J., Moreno-Arenas, G., 2019.  On residuals in generalized johnson sb regressions.  Applied Mathematical Modelling 67, 62–73.
  * Lemonte and Patriota (2011) Lemonte, A.J., Patriota, A.G., 2011.  Multivariate elliptical models with general parameterization.  Statistical Methodology 8, 389–400.
  * Li et al. (2022) Li, D., Al-Mahamda, M.F., Song, Y., Feng, S., Sze, N.N., 2022.  An alternate crash severity multicategory modeling approach with asymmetric property.  Analytic Methods in Accident Research 35, 100218.
  * Maier (2014) Maier, M.J., 2014.  Dirichletreg: Dirichlet regression for compositional data in r.  Research Report Series / Department of Statistics and Mathematics 125.
  * Maski et al. (2021) Maski, K.P., Colclasure, A., Little, E., Steinhart, E., Scammell, T.E., Navidi, W., Diniz Behn, C., 2021.  Stability of nocturnal wake and sleep stages defines central nervous system disorders of hypersomnolence.  Sleep 44, zsab021.
  * Melo et al. (2022) Melo, T.F., Vargas, T.M., Lemonte, A.J., Moreno-Arenas, G., 2022.  Higher-order asymptotic refinements in the multivariate dirichlet regression model.  Communications in Statistics-Simulation and Computation 51, 53–71.
  * Morais et al. (2018) Morais, J., Thomas-Agnan, C., Simioni, M., 2018.  Using compositional and dirichlet models for market share regression.  Journal of Applied Statistics 45, 1670–1689.
  * Ohayon et al. (2017) Ohayon, M., Wickwire, E.M., Hirshkowitz, M., Albert, S.M., Avidan, A., Daly, F.J., Dauvilliers, Y., Ferri, R., Fung, C., Gozal, D., et al., 2017.  National sleep foundation’s sleep quality recommendations: first report.  Sleep health 3, 6–19.
  * Pereira et al. (2020) Pereira, G.H., Scudilio, J., Santos-Neto, M., Botter, D.A., Sandoval, M.C., 2020.  A class of residuals for outlier identification in zero adjusted regression models.  Journal of Applied Statistics 47, 1833–1847.
  * Pereira (2019) Pereira, G.H.A., 2019.  On quantile residuals in beta regression.  Communications in Statistics-Simulation and Computation 48, 302–316.
  * Reichert et al. (2022) Reichert, S., Berger, V., Dos Santos, D.J.F., Lahdenperä, M., Nyein, U.K., Htut, W., Lummaa, V., 2022.  Age related variation of health markers in asian elephants.  Experimental Gerontology 157, 111629.
  * Ritmala-Castren et al. (2015) Ritmala-Castren, M., Virtanen, I., Leivo, S., Kaukonen, K.M., Leino-Kilpi, H., 2015.  Sleep and nursing care activities in an intensive care unit.  Nursing & health sciences 17, 354–361.
  * Ross (2020) Ross, G.J., 2020.  Tracking the evolution of literary style via dirichlet–multinomial change point regression.  Journal of the Royal Statistical Society: Series A (Statistics in Society) 183, 149–167.
  * Schaid et al. (2019) Schaid, D.J., Tong, X., Batzler, A., Sinnwell, J.P., Qing, J., Biernacka, J.M., 2019.  Multivariate generalized linear model for genetic pleiotropy.  Biostatistics 20, 111–128.
  * Scudilio and Pereira (2020) Scudilio, J., Pereira, G.H., 2020.  Adjusted quantile residual for generalized linear models.  Computational Statistics 35, 399–421.
  * Sen et al. (2010) Sen, P.K., Singer, J.M., de Lima, A.C.P., 2010.  From finite sample to asymptotic methods in statistics. volume 28.  Cambridge University Press.
  * Stephens (1986) Stephens, M.A., 1986.  Tests based on edf statistic, in: Goodness-of-fit Techniques. Marcel Dekker, pp. 97–193.
  * Tempesta et al. (2018) Tempesta, D., Socci, V., De Gennaro, L., Ferrara, M., 2018.  Sleep and emotional processing.  Sleep medicine reviews 40, 183–195.
  * Tsagris and Stewart (2018) Tsagris, M., Stewart, C., 2018.  A dirichlet regression model for compositional data with zeros.  Lobachevskii Journal of Mathematics 39, 398–412.
  * Tsilimigras and Fodor (2016) Tsilimigras, M.C., Fodor, A.A., 2016.  Compositional data analysis of the microbiome: fundamentals, tools, and challenges.  Annals of epidemiology 26, 330–335.
  * Veje et al. (2021) Veje, M., Studahl, M., Thunström, E., Stentoft, E., Nolskog, P., Celik, Y., Peker, Y., 2021.  Sleep architecture, obstructive sleep apnea and functional outcomes in adults with a history of tick-borne encephalitis.  PLoS One 16, e0246767.
  * Vieira and Salmon (2019) Vieira, B.H., Salmon, C.E.G., 2019.  A principled multivariate intersubject analysis of generalized partial directed coherence with dirichlet regression: Application to healthy aging in areas exhibiting cortical thinning.  Journal of Neuroscience Methods 311, 243–252.
  * Yoo et al. (2022) Yoo, J., Sun, Z., Ma, Q., Chung, D., Kim, Y.M., 2022.  A guideline for the statistical analysis of compositional data in immunology.  arXiv preprint arXiv:2201.07945 .
  * Zhao et al. (2019) Zhao, D., Westfall, J., Coulston, J.W., Lynch, T.B., Bullock, B.P., Montes, C.R., 2019.  Additive biomass equations for slash pine trees: Comparing three modeling approaches.  Canadian Journal of Forest Research 49, 27–40.

## Appendix

In this Appendix, we prove Theorem 2.2.1 using the definitions presented in
Section 2. Under the correct model, if $\tau_{i}$ is consistently estimated,
$y_{i}^{\\{b\\}}$, for $b=1,2,\ldots,B$, and $y_{i}$ have the same asymptotic
distribution. As a result $l_{i}^{\\{b\\}}$ and $l_{i}$ have also the same
asymptotic distribution and

$a_{i}\xrightarrow{D}U_{D}(0;B),$ (5)

where $U_{D}(0;B)$ denotes the discrete uniform distribution with parameters 0
and $B$.

If $0<k<1/(B+1)$, then

$P(u_{i}\leq k)=P(u_{i}\leq k,a_{i}=0)=P(u_{i}\leq
k|a_{i}=0)P(a_{i}=0)=P(U^{\\{1\\}}\leq k)P(a_{i}=0),$

where $U^{\\{1\\}}$ is a uniform distributed random variable with parameters 0
and $1/(B+1)$. Therefore, considering (5),

$P(u_{i}\leq k)\xrightarrow{P}\frac{k}{1/(B+1)}\,\frac{1}{B+1}=k.$ (6)

If $1/(B+1)\leq k<1$, then

$P(u_{i}\leq
k)=\sum_{j=0}^{\lambda-1}P(a_{i}=j)+P\left(\frac{\lambda}{B+1}<u_{i}\leq
k,a_{i}=\lambda\right),$

where $\lambda$ is the largest integer less or equal to $k(B+1)$. Then,

$\begin{split}P(u_{i}\leq
k)&=\sum_{j=0}^{\lambda-1}P(a_{i}=j)+P\left(\frac{\lambda}{B+1}<u_{i}\leq
k|a_{i}=\lambda\right)P(a_{i}=\lambda)\\\
&=\sum_{j=0}^{\lambda-1}P(a_{i}=j)+P(U^{\\{2\\}}\leq
k)P(a_{i}=\lambda),\end{split}$

where $U^{\\{2\\}}$ is a uniform distributed random variable with parameters
$\lambda/(B+1)$ and $(\lambda+1)/(B+1)$. Therefore, considering (5),

$P(u_{i}\leq
k)\xrightarrow{P}\frac{\lambda}{B+1}+\frac{k-\lambda/(B+1)}{1/(B+1)}\,\frac{1}{B+1}=\frac{\lambda}{B+1}+\frac{k(B+1)-\lambda}{B+1}=k.$
(7)

From (6) and (7), $P(u_{i}\leq k)\xrightarrow{P}k$, for all $0<k<1$ and so

$u_{i}\xrightarrow{D}U(0;1).$ (8)

Using (8) and the Sverdrup’s theorem (Sen et al., 2010, page 180),

$\Phi^{-1}\\{u_{i}\\}\xrightarrow{D}\Phi^{-1}\\{U(0;1)\\},$

and so

$s_{i}\xrightarrow{D}N(0;1).$


11institutetext: PES University, Bangaluru KA 560085, India 22institutetext:
Regional Remote Sensing Centre - South, ISRO, Bengaluru KA 560037, India
22email<EMAIL_ADDRESS>
# Analysis and application of multispectral data for water segmentation using
machine learning

Shubham Gupta Work done partially while interning at CDSAML, PES University
and RRSC-S, ISRO.11 Uma D 11 Hebbar  R 22

###### Abstract

Monitoring water is a complex task due to its dynamic nature, added
pollutants, and land build-up. The availability of high-resolu-tion data by
Sentinel-2 multispectral products makes implementing remote sensing
applications feasible. However, overutilizing or underutilizing multispectral
bands of the product can lead to inferior performance. In this work, we
compare the performances of ten out of the thirteen bands available in a
Sentinel-2 product for water segmentation using eight machine learning
algorithms. We find that the shortwave infrared bands (B11 and B12) are the
most superior for segmenting water bodies. B11 achieves an overall accuracy of
$71\%$ while B12 achieves $69\%$ across all algorithms on the test site. We
also find that the Support Vector Machine (SVM) algorithm is the most
favourable for single-band water segmentation. The SVM achieves an overall
accuracy of $69\%$ across the tested bands over the given test site. Finally,
to demonstrate the effectiveness of choosing the right amount of data, we use
only B11 reflectance data to train an artificial neural network, BandNet. Even
with a basic architecture, BandNet is proportionate to known architectures for
semantic and water segmentation, achieving a $92.47$ mIOU on the test site.
BandNet requires only a fraction of the time and resources to train and run
inference, making it suitable to be deployed on web applications to run and
monitor water bodies in localized regions. Our codebase is available at
https://github.com/IamShubhamGupto/BandNet.

###### Keywords:

Water Sentinel-2 machine-learning artificial neural network BandNet

##  1 Introduction

Water is one of the most essential resources on our planet. With the rise in
population in cities, water bodies have shrunk or disappeared entirely. In the
coastal regions, due to global warming, the increase in water levels has begun
to submerge land in low-lying areas[13]. To take action and prevent
irreparable damages, spatio-temporal awareness of water is a necessity.

The recent advancements in Convolutional Neural Network (CNN) models show
promising results in object classification[22, 14] and segmentation[5, 12, 20]
from images and reflectance values[21, 15, 11]. The reflectance models rely on
the combination of either or all of the visible spectrum, near-infrared, and
shortwave-infrared spectrum bands. To our knowledge, there is no ranking of
the multispectral bands for their usability for water segmentation. In this
work, we study the relationship between water segmentation performance and
multispectral bands using several machine learning algorithms[9, 2, 23, 17, 1,
6, 18, 8]. The goal is to find the multi-spectral band and machine learning
algorithm that is best suited to segmenting water bodies.

By utilizing the best multispectral band to delineate water bodies, we show
that a single band is adequate to perform similar to large segmentation
models[5] trained on multispectral false-infrared images  that are composed of
B8, B4, and B3. We achieve this result on a simple Artificial Neural Network,
BandNet. Unlike existing literature, BandNet does not treat multispectral data
as images but rather sees each pixel as a unique data-point. Furthermore, it
requires a fraction of the training time and resources to produce respectable
results. To summarize our contributions to this paper:

  1. 1.

We compare the water segmentation performance of the multispectral bands
present in the Sentinel-2 product using multiple machine learning algorithms.

  2. 2.

We compare the performances of various machine learning algorithms for single-
band water segmentation.

  3. 3.

We demonstrate the importance of data selection by training a simple ANN,
BandNet, that performs similar to other segmentation models - trained on
images and relfectance data.

##  2 Satellite Data and Study Site

This study uses the multi-spectral products from the European Space Agency
satellite constellation Sentinel-2. All the products are made available in L2A
mode, giving us direct access to Bottom of Atmosphere (BOA) reflectance[3]. We
study bands B2, B3, B4, B5, B6, B7, B8, B8A, B11, and B12 which are re-sampled
to $10m$ spatial resolution using SNAP, SeNtinel Applications Platform.

Since one of the requirements of training Deeplabv3+[5] is an extensively
annotated dataset, we required images of water bodies from more than one
Sentinel-2 product. For this reason, we acquired a total of eight Sentinel-2
products. As for the image from the products themselves, we utilize false
infrared images generated from B8, B4, and B3 of each product. The entire
image is 10980x10980 pixels. Since we cannot process an image of this scale at
once, we break down the false-infrared images into 549x549 pixel tiles.

Table 1: Metadata of the Sentinel-2 data used. *Rotation is applied to images as part of data augmentation. Tile |  Date |  Type |  Images
---|---|---|---
T43PGQ |  2019-01-04 |  Lakes |  186*
T44RQN |  2019-01-21 |  Rivers |  64
T43PFS |  2019-02-16 |  Lakes, Rivers |  47
T43QGU |  2019-03-18 |  Lakes, Rivers |  172*
T43PFQ |  2019-03-28 |  Lakes |  310*
T44QLE |  2019-03-30 |  Rivers |  175
T43RFJ |  2019-04-20 |  Rivers |  242*
T43QHV |  2019-04-24 |  Lakes |  250*

We use the scene classification map available within each Sentinel-2 product
to annotate water bodies. We discard tiles and the corresponding annotation if
they do not contain any water bodies. The resulting data acquired from these
products is summarized in Table 1. We obtain 1446 images with annotations
after applying a simple rotation augmentation to five of the eight products.

(a)

(b)

(c)

(d)

Figure 1: False infrared image generated from the Sentinel-2 product T43PGQ
subset dated 2019-01-04. (a)Preview of subset1. (b), (d) Corresponding water
body annotation generated from scene classification map. (c)Preview of
subset2.

On the other hand, a single Sentinel-2 product can provide $120$ million data
points as the values are based on each pixel. Due to hardware limitations, we
can only take a subset of the product at a time. Figure1 shows the false-
infrared images of the two subsets we will be using to evaluate our findings.
Image (a) in Figure1 represents a subset with North latitude 12.963, West
longitude 77.630, South latitude 12.889, East longitude 77.703. This subset
contains Bellandur lake and has a good representation of land built-up, water
bodies, and vegetation. Image (c) in Figure1 represents a subset with North
latitude 12.655, West longitude 77.131, South latitude 12.581, East longitude
77.205. This subset contains vegetated land and lakes at close proximity. We
will be referring to them as subset1 and subset2 respectively.

##  3 Methodology Used

###  3.1 Data processing

From subset1, we consider an equal number of water and non-water data points.
For reproducibility, we set a seed value[16] so that the same data points will
be selected when the experiment is re-run. The dataset is split into test,
train, and validation sets, each being mutually exclusive. We evaluate and
report the performance of the machine learning algorithms on the validation
set. The validation set has never been seen before by the algorithms. Data
points from subset2 are entirely used to evaluate the algorithms and models.

###  3.2 Band reflectance analysis

We employ eight statistical machine learning algorithms on the reflectance
data of individual bands to gauge their performance for water segmentation,
namely: Logistic Regression (LR)[9], Gaussian Naive Bayes (GNB)[2],
K-neighbours (KN)[23], Decision Tree (DT)[17], Random Forest (RF)[1], XGBoost
(XGB)[6], Stochastic Gradient Descent (SGD)[18], and Support Vector Machine
(SVM)[8]. All eight algorithms are first made to fit the training set,
followed by testing their performance on the validation set.

###  3.3 BandNet

Based on the performance of individual bands, further discussed in Section5,
we designed a simple ANN architecture, BandNet, to create segmentation maps of
water bodies using raw reflectance data as input.

Figure 2: Proposed architecture of BandNet for segmentation of Water Bodies.

BandNet consists of two Dense layers with ReLU activation functions. Before
the classification head, a Dropout layer of 0.4 is in place to deal with the
class imbalance problem and improve generalization. The classification head is
a Dense layer with Sigmoid activation. The architecture of BandNet is shown in
Figure2.

###  3.4 Multi-spectral Image Analysis

We generate false colour infrared images from B8, B4, and B3 to train
Deeplabv3+[5], a deep neural network architecture that uses Atrous Spatial
Pyramid Pooling (ASPP)[4] and a encoder–decoder[12] structure. The ASPP is
used to encode the contextual information with the help of varying atrous
rates $r$.

$y[i]=\sum_{k}x[i+r.k]w[k]$ (1)

The decoder then uses the contextual information to generate fine object
boundaries. Finally, depthwise separable convolution[7] is applied to both
ASPP and the decoder module for reducing computational complexity while
preserving the performance of the model.

##  4 Implementation details

Here we define the details used to train the various models and algorithms. We
use an Nvidia GTX 1050 graphics card, 16GB of memory, and an Intel i5 8300H
processor for training and inference. The SVM classifier uses a radial basis
function (rbf) kernel. The RF classifier uses 100 estimators. The KN
classifier uses seven neighbors to vote. The SGD classifier uses a Huber loss
function with an L1 penalty for 25 iterations. The remaining machine learning
algorithms ran with no extra parameters.

For each multispectral band and model, we generate a confusion matrix, which
then helps us calculate the mean intersection over union (mIoU). mIoU for N
classes can be defined in terms of true positives (TP) as:

$IoU\textsubscript{i}=\frac{TP\textsubscript{i}}{SumOfRows\textsubscript{i}+SumOfColumns\textsubscript{i}-TP\textsubscript{i}}$
(2) $mIoU=\frac{\sum_{n=1}^{N}IoU\textsubscript{i}}{N}$ (3)

We train BandNet with a BinaryCrossEntropy loss function and BinaryAccuracy as
the metric. The optimizer Adam has a learning rate of $1e-4$. We use
EarlyStopping to avoid overfitting with a patience value of 5.

We run the DeepwatermapV2[11] and WatNet[15] directly using pretrained weights
for inference on B2, B3, B4, B8A, B11, and B12 based on the guides provided by
the authors.

The Deeplabv3+ uses a learning rate of $0.001$, weight decay of $0.00004$,
atrous rates at (6,12,18), output stride of $16$, momentum of $0.9$, decoder
output stride at $4$, batch size of $1$ for 30K iterations. We use
Xception65[7] and MobileNetv2[19] as the network backbones pretrained on the
PASCALVOC-2012 dataset[10].

##  5 Results and Discussion

Table 2: Comparison of single band water segmentation performance (mIoU) using various machine learning models on validation set of subset1. The mIoU is colour coded from low to high as brick-red to blue. The Percent row and column describes the accuracy of a band across various algorithms and the accuracy of an algorithm across various bands respectively. The colour coding from low to high follows as  red to  teal. Band |  LR |  GNB |  RF |  KN |  DT |  SGD |  XGB |  SVM |  |  Percent
---|---|---|---|---|---|---|---|---|---|---
B02 |  [HTML]a4180b36.93 |  [HTML]9e0b0534.23 |  [HTML]b1321842.07 |  [HTML]b2361942.74 |  [HTML]b2341942.45 |  [HTML]a0100735.26 |  [HTML]b2351942.62 |  [HTML]ae2c1540.76 |  |  [HTML]e67c730.35
B03 |  [HTML]a4180b36.79 |  [HTML]99000032.04 |  [HTML]b0301741.58 |  [HTML]b0311741.69 |  [HTML]b0321741.89 |  [HTML]ae2d1540.99 |  [HTML]af2f1641.45 |  [HTML]b1321842.06 |  |  [HTML]e67c730.35
B04 |  [HTML]b63e1d44.37 |  [HTML]b2351942.66 |  [HTML]c7622e51.34 |  [HTML]c8632f51.58 |  [HTML]c8642f51.82 |  [HTML]c0522748.33 |  [HTML]c7622e51.43 |  [HTML]df934661.06 |  |  [HTML]ec9d970.43
B05 |  [HTML]c2572949.17 |  [HTML]c0522748.23 |  [HTML]d0743755.02 |  [HTML]d1763855.3 |  [HTML]d2793955.86 |  [HTML]d8843f58.11 |  [HTML]d47c3b56.51 |  [HTML]d47c3b56.48 |  |  [HTML]f0b0ab0.47
B06 |  [HTML]facc6172.26 |  [HTML]c7b87077.27 |  [HTML]dec46c76.02 |  [HTML]fbce6272.68 |  [HTML]d1bd6e76.75 |  [HTML]a0a37779.44 |  [HTML]a1a47779.41 |  [HTML]aaa87678.92 |  |  [HTML]eff9f40.67
B07 |  [HTML]fdd36473.66 |  [HTML]9fa27879.53 |  [HTML]d5c06e76.52 |  [HTML]fed46573.88 |  [HTML]c9b97077.16 |  [HTML]b1ac7478.5 |  [HTML]c1b57177.62 |  [HTML]b7af7378.18 |  |  [HTML]e3f4eb0.67
B08 |  [HTML]f5c15c70.15 |  [HTML]efcd6975.08 |  [HTML]fbd46674.4 |  [HTML]bab17378.02 |  [HTML]e3c76b75.72 |  [HTML]fcd16373.16 |  [HTML]a7a77679.05 |  [HTML]9aa07979.78 |  |  [HTML]fefcfc0.66
B8A |  [HTML]e1c66b75.81 |  [HTML]64838382.77 |  [HTML]c9b97077.16 |  [HTML]fbcf6272.82 |  [HTML]bfb47277.72 |  [HTML]b7af7378.19 |  [HTML]bab17378.03 |  [HTML]9ca17879.69 |  |  [HTML]ccebdb0.68
B11 |  [HTML]89977c80.75 |  [HTML]0b539487.71 |  [HTML]778d7f81.71 |  [HTML]afab7578.61 |  [HTML]6d888182.3 |  [HTML]bab17378 |  [HTML]6f898182.2 |  [HTML]43718984.65 |  |  [HTML]57bb8a0.71
B12 |  [HTML]f8c75f71.33 |  [HTML]26628e86.21 |  [HTML]bab17378.01 |  [HTML]b6af7378.21 |  [HTML]959d7a80.06 |  [HTML]c6b77077.36 |  [HTML]42718984.68 |  [HTML]42718984.69 |  |  [HTML]96d5b60.69
Percent |  [HTML]e67c730.61 |  [HTML]fbeeed0.65 |  [HTML]f0f9f50.65 |  [HTML]fbeeed0.65 |  [HTML]cbeadb0.66 |  [HTML]fcf4f30.65 |  [HTML]96d5b60.67 |  [HTML]57bb8a0.69 |  |

In table2, we present our results of generating water segmentation from each
of the ten bands paired with each of the eight machine learning algorithms. We
calculate the percent column for a band as a fraction of the sum of the
achieved mIoU across all algorithms, to the maximum attainable mIoU. The
percent row for each algorithm is calculated as the sum of the achieved mIoU
across all bands divided by the maximum attainable mIoU.

In that order, B11, B12, and B8A are the best-performing bands across all
tested algorithms. B11 outperforms B2 by an absolute value of $0.36$ which
makes B11 $102\%$ better than B2 for water segmentation. We also notice a
trend of SWIR bands (B11 and B12), followed by NIR bands (B8, B8A) and finally
visible spectrum bands (B2, B3, B4) in terms of water segmentation
performance. However, the algorithms used to test these bands do not present
such large variations in performance. The best performing algorithm, Support
Vector Machines[8], outperforms Linear Regression[9] by an absolute value of
$0.08$ or $13.11\%$.

Table 3: Comparison of BandNet to other models on subset1 in terms of performance, training time, parameters. mIOU: mean Intersection over Union. Time is measured in Hours. $\downarrow$: The lower the better. $\uparrow$: The higher the better.  Model |  Data |  Parameters $\downarrow$ |  Time $\downarrow$ |  mIOU $\uparrow$
---|---|---|---|---
DeepWaterMapv2[11] |  B2-B3-B4-B8A-B11-B12 |  37.2M |  - |  89.01
WatNet[15] |  B2-B3-B4-B8A-B11-B12 |  3.4M |  - |  89.74
BandNet (Ours) |  B8-B4-B3 |  1217 |  0.41 |  89.92
MobileNetv2[5, 19] |  Image |  2.1M |  12 |  90.27
Xception65[5, 7] |  Image |  41M |  16 |  92.44
BandNet (Ours) |  B11 |  1153 |  0.45 |  92.47
BandNet (Ours) |  B12-B11 |  1185 |  0.27 |  92.96
BandNet (Ours) |  B12-B11-B8A |  1217 |  0.65 |  92.99
Table 4: Comparison of BandNet to reflectance models on subset1 in terms of performance, training time, parameters. mIOU: mean Intersection over Union. Time is measured in Hours. $\uparrow$: The higher the better. Model |  Data |  mIOU $\uparrow$
---|---|---
WatNet[15] |  B2-B3-B4-B8A-B11-B12 |  86.92
BandNet (Ours) |  B12-B11-B8A |  87.00
BandNet (Ours) |  B12-B11 |  87.33
BandNet (Ours) |  B11 |  87.42
DeepWaterMapv2[11] |  B2-B3-B4-B8A-B11-B12 |  89.33

In table 3, we compare the performance of BandNet on different band
combinations to existing solutions for water segmentation.This comparison is
carried out on subset1. We directly use Deepwatermapv2 and WatNet using their
existing weights provided by the authors. Both models perform respectably
considering they have never seen the data before in training or testing
phases.

To make a fair comparison of images to reflectance data, we compared Deeplab
V3+ to BandNet trained on the reflectance of B8-B4-B3. We notice BandNet
performs similarly to Deeplabv3+ with MobileNetv2 as the backbone. However,
BandNet trains at a fraction of the parameters and time required by Deeplabv3+
for the same performance. This performance gap is bridged when we train
BandNet on B11, the best performing multispectral band. We observe that
BandNet with B11 outperforms Deeplabv3+ by an absolute value of $2.2$ and
$0.03$ mIOU on MobileNetv2 and Xception65 backbones, respectively. BandNet
further surpasses Deeplab3+ if we pass along additional high-performing bands
with B11.

In table4, we compare the inference of models on subset2. BandNet sees a
slight drop in performance but is still comparable to Deepwatermapv2 and
WatNet. These results underscore the idea behind this work that choosing the
right set of data is important when building predictive models. We believe
this drop in performance is due to the simplistic nature of BandNet and that
it was originally designed for one multispectral band only. However, this
makes it convenient to re-train and deploy BandNet over a very localized
region using minimal resources and time. The small size of BandNet allows it
to be hosted on a web application and run inferences at a very low cost. This
is extremely helpful if someone were to host their own water body monitoring
solution.

(a)

(b)

(c)

Figure 3: Application of BandNet: Monitoring of Bellandur lake using BandNet
in subset1. The maps are generated on the T43PGQ subset dated: (a) 2019-01-04,
(b) 2020-03-29 and (c) 2021-03-04.

##  6 Conclusion

In this paper, we have compared the performances of individual bands for water
segmentation, using machine learning algorithms. We create a hierarchy of
multispectral bands based on their performance, placing SWIR bands (B11, B12)
on the top, followed by NIR bands (B8, B8A) and finally visible spectrum bands
(B2, B3, B4). We also observed that the Support Vector Machine followed by XG
Boost algorithm are favourable is single band water segmentation compared to
other algorithms.

Using the best performing band, B11, we developed a simple ANN that is able to
compete with specialized segmentation architectures in performance while
requiring only a fraction of the time and resources to train. This lightweight
nature of BandNet makes it suitable for deploying it on web applications to
monitor water bodies in localized regions.

The objective we want to underscore in this study is that using the right
amount and type of data is sufficient to compete with existing solutions with
no improvements to the architecture. Now that we have established a hierarchy
of multispectral bands, we can further look into designing a novel
architecture or training regime for water segmentation. Our process can also
be extended to other features such as vegetation or built-up land and can be
used in combinations for multi-class segmentation.

####  6.0.1 Acknowledgment

This work has been supported by Center of Data Science and Applied Machine
Learning, Computer Science and Engineering Department of PES University, and
Regional Remote Sensing Centre - south. We would like to thank Dr. Shylaja,
S.S. of PES University and Dr. K. Ganesha Raj of Regional Remote Sensing
Centre - south for the opportunity to carry out this work.

## References

  * [1] Breiman, L.: Random forests. Machine Learning 45(1), 5–32 (2001). https://doi.org/10.1023/A:1010933404324, http://dx.doi.org/10.1023/A%3A1010933404324
  * [2] Chan, Golub, L.: Updating formulae and a pairwise algorithm for computing sample variances. http://i.stanford.edu/pub/cstr/reports/cs/tr/79/773/CS-TR-79-773.pdf (1979)
  * [3] Chen, J., Li, Y., Ma, Q., Shen, X., Zhao, A., Li, J.: Preliminary evaluation of sentinel-2 bottom of atmosphere reflectance using the 6sv code in beijing area. In: IGARSS 2018 - 2018 IEEE International Geoscience and Remote Sensing Symposium. pp. 7760–7763 (2018). https://doi.org/10.1109/IGARSS.2018.8517598
  * [4] Chen, L.C., Papandreou, G., Schroff, F., Adam, H.: Rethinking atrous convolution for semantic image segmentation (06 2017)
  * [5] Chen, L.C., Zhu, Y., Papandreou, G., Schroff, F., Adam, H.: Encoder-decoder with atrous separable convolution for semantic image segmentation (02 2018)
  * [6] Chen, T., Guestrin, C.: XGBoost: A scalable tree boosting system. In: Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. pp. 785–794. KDD ’16, ACM, New York, NY, USA (2016). https://doi.org/10.1145/2939672.2939785, http://doi.acm.org/10.1145/2939672.2939785
  * [7] Chollet, F.: Xception: Deep learning with depthwise separable convolutions. In: 2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR). pp. 1800–1807 (2017)
  * [8] Cortes, C., Vapnik, V.: Support-vector networks. Machine learning 20(3), 273–297 (1995)
  * [9] Cox, D.R.: The regression analysis of binary sequences. Journal of the Royal Statistical Society: Series B (Methodological) 20(2), 215–232 (1958)
  * [10] Everingham, M., Van Gool, L., Williams, C.K.I., Winn, J., Zisserman, A.: The PASCAL Visual Object Classes Challenge 2012 (VOC2012) Results. http://www.pascal-network.org/challenges/VOC/voc2012/workshop/index.html
  * [11] Isikdogan, L., Bovik, A., Passalacqua, P.: Seeing through the clouds with deepwatermap. IEEE Geoscience and Remote Sensing Letters PP, 1–5 (11 2019). https://doi.org/10.1109/LGRS.2019.2953261
  * [12] Islam, M., Rochan, M., Bruce, N., Wang, Y.: Gated feedback refinement network for dense image labeling. pp. 4877–4885 (07 2017). https://doi.org/10.1109/CVPR.2017.518
  * [13] Kazama, S., Oki, T.: The effects of climate change on water resources (09 2006)
  * [14] Liu, Z., Hu, H., Lin, Y., Yao, Z., Xie, Z., Wei, Y., Ning, J., Cao, Y., Zhang, Z., Dong, L., et al.: Swin Transformer V2: Scaling Up Capacity and Resolution. arXiv preprint arXiv:2111.09883 (2021)
  * [15] Luo, X., Tong, X., Hu, Z.: An applicable and automatic method for earth surface water mapping based on multispectral images. International Journal of Applied Earth Observation and Geoinformation 103, 102472 (2021)
  * [16] Picard, D.: Torch. manual_seed (3407) is all you need: On the influence of random seeds in deep learning architectures for computer vision. arXiv preprint arXiv:2109.08203 (2021)
  * [17] Quinlan, J.R.: Induction of decision trees. Machine Learning 1, 81–106 (1986)
  * [18] Ruder, S.: An overview of gradient descent optimization algorithms. arXiv preprint arXiv:1609.04747 (2016)
  * [19] Sandler, M., Howard, A., Zhu, M., Zhmoginov, A., Chen, L.C.: Mobilenetv2: Inverted residuals and linear bottlenecks. pp. 4510–4520 (06 2018). https://doi.org/10.1109/CVPR.2018.00474
  * [20] Wei, Y., Hu, H., Xie, Z., Zhang, Z., Cao, Y., Bao, J., Chen, D., Guo, B.: Contrastive learning rivals masked image modeling in fine-tuning via feature distillation. arXiv preprint arXiv:2205.14141 (2022)
  * [21] Yang, L., Driscol, J., Sarigai, S., Wu, Q., Lippitt, C.D., Morgan, M.: Towards synoptic water monitoring systems: A review of ai methods for automating water body detection and water quality monitoring using remote sensing. Sensors 22(6) (2022). https://doi.org/10.3390/s22062416, https://www.mdpi.com/1424-8220/22/6/2416
  * [22] Zhang, H., Li, F., Liu, S., Zhang, L., Su, H., Zhu, J., Ni, L.M., Shum, H.Y.: Dino: Detr with improved denoising anchor boxes for end-to-end object detection. arXiv preprint arXiv:2203.03605 (2022)
  * [23] Zhang, M.L., Zhou, Z.H.: A k-nearest neighbor based algorithm for multi-label classification. In: Hu, X., Liu, Q., Skowron, A., Lin, T.Y., Yager, R.R., Zhang, B. (eds.) GrC. pp. 718–721. IEEE (2005), http://cs.nju.edu.cn/zhouzh/zhouzh.files/publication/grc05.pdf