# A new approach to the gender pay gap decomposition by economic
activity111Supported by the MINECO grants MTM2017-82724-R, MTM2015-71217-R,
and by the Xunta de Galicia (Grupos de Referencia Competitiva ED431C-2016-015
and Centro Singular de Investigación de Galicia ED431G/01), all of them
through the ERDF.

María José Lombardía1 Esther López-Vizcaíno2
Cristina Rueda3
1 Universidade da Coruña, CITIC, Spain,
2 Instituto Galego de Estatística, Spain,
3Universidad de Valladolid, Spain

Abstract

The aim of this paper is to present an original approach to estimate the
gender pay gap. We propose a model-based decomposition, similar to the most
popular approaches, where the first component measures differences in group
characteristics and the second component measures the unexplained effect; the
latter being the real gap. The novel approach incorporates model selection and
bias correction. The pay gap problem in a small area context is considered in
this paper, although the approach is flexible to be applied to other contexts.

Specifically, the methodology is validated for analysing wage differentials by
economic activities in the region of Galicia (Spain) and by analysing
simulated data from an experimental design that imitates the generation of
real data. The good performance of the proposed estimators is shown in both
cases, specifically when compared with those obtained from the widely used
Oaxaca-Blinder approach.

##  1 Introduction

Gender wage inequality is a documented fact that happens in almost all of the
industrialized countries (Juhn et al., 1993, Oaxaca and Ransom, 1994,1998,
Aláez and Ullibarri, 2001, Moral-Arce et al., 2011, Fortin et al., 2011) and
interest in the quantification of this gap has increased in recent years.

In the case of the Spanish labor market, a number of actions were taken to
provide women with real access to employment with full social and economic
rights, with special emphasis on the reconciliation of family and working
life. This and other reasons led a massive incorporation of women into the
labor market, but it has not resulted in equality, that is, women’s
unemployment rates are higher in most sectors, the jobs that women take do not
involve the same degree of responsibility or decision-making power and women’s
participation is limited to a few sectors of the economy. But more important
is that, for a similar job, it has been shown that men have better wages than
women (Fernández et al., 2000; Aláez and Ullibarri, 2001; Aláez et al.,
2000-2003; De la Rica et al., 2008; Moral-Arce et al., 2011).

The gender pay gap ($GPG$), calculated as the difference between men’s average
hourly earnings and women’s average hourly earnings, usually expressed as a
percentage of the men’s average hourly earnings, is widely used as the key
indicator to study the progress in the European Union (EU), see Leythienne and
Ronkowski (2018), among others. Moreover, most researchers in this field
(Jann, 2008, Moral-Arce et al., 2011, Anastasiade and Tillé, 2017a, 2017b,
Hlavac, 2018) work with the transformed variable to logarithms, $Y=log(W)$,
where $W$ is the hourly wages. The rationale behind the use of the logarithmic
transformation is that the salaries have a skew distribution. However, the
skewness is also exhibited on the log scale. Finding the optimal
transformation is not straightforward. As an alternative, Graf et al. (2019)
propose considering a much more flexible distribution called generalized beta
of the second kind (GB2). In this work, we follow the ideas of most
researchers. We are aware that it would be interesting to incorporate this
proposal, but it goes beyond the objectives of this paper.

Typically the sample estimators are used directly as the GPG measure because
the sample size are often very high. Instead, a more formal definition of the
GPG is used here, as we deal with scenarios with small sample sizes. It is
given, using expectations, as follows:

$GPG=\frac{E(W_{m})-E(W_{w})}{E(W_{m})},$

where the subscripts $m$ are the men’s wage and $w$ the women’s wage.

The most widely used approach to analyze the gender gap is the decomposition
due to Oaxaca (1973) and Blinder (1973); Specifically, Eurostat applies this
approach to decompose the estimator of $GPG$ (Leythienne and Ronkowski, 2018).
The Oaxaca-Blinder (OB) decomposition breaks down the difference
$\Delta=E(Y_{m})-E(Y_{w})$ into two components: $\Delta=Q+U$, where the first
component explains the difference between the observed productive
characteristics, such as education or work experience, and the second accounts
for differences in the structure of the model (Anastasiade and Tillé, 2017a,
2017b). The U component is known as the unexplained component and is usually
considered to be the wage discrimination by gender in the labor market. By
instances, Popli (2013) states that the unexplained wage gap includes the
effect of labour market discrimination, unobservable variables and omitted
variables. The methodology presented in this paper accounts for the bias from
unobservable or omitted variables.

So, the $GPG$ decomposition is:

${GPG}_{Q}={GPG}\frac{{Q}}{\Delta}\ \quad\ \mbox{and}\quad\
{GPG}_{U}={GPG}\frac{{U}}{\Delta},$ (1)

which represent the explained and unexplained part of the $GPG$, respectively.
The OB method, including the definitions of $Q$ and $U$, is briefly revised in
Section 2.

Despite its widespread use, the standard Oaxaca-Blinder decomposition has
important limitations; firstly, too simple parametric functions are often used
for the regression model; secondly, the results may depend heavily on the set
of explanatory variables selected; and thirdly, the unexplained variability is
not always taken into account and generates bias.

This paper aims to adapt the methodology to estimate $GPG_{Q}$ and $GPG_{U}$
in small areas, which are disaggregated levels of the population, and
specifically analyze the wage differentials by economic activities in the
region of Galicia (Spain). Often in small area contexts, the sample sizes are
too small, which implies the standard estimators are unreliable. In
particular, in our application for 25% of the economic activities, we have
information on less than 26 male workers, being the sample sizes for women
even smaller.

In fact, no official estimators of wages are provided at the economic activity
level because direct estimates have very low accuracy. The latter is a typical
SAE (Small Area Estimation) problem. The most extended SAE approach is to give
a model-based estimator, where auxiliary information in terms of explanatory
variables is incorporated. Besides, the model is often defined as a mixed
effect model where a random effect represents the areas (see Rao and
Molina(2015).

The proposal here is to combine SAE methodology and OB methodology by deriving
a decomposition inspired by (1), adapted to a small area scenario in such a
way that for a given area $d$, $GPG_{d}$ is also decomposing into the
$GPG_{Qd}$ and $GPG_{Ud}$. These components are efficiently estimated using an
approach in three steps that use SAE techniques and overcome OB’s drawback,
commented above. First, by doing model selection using an AIC criterion, and
second, by including a bias correction term that accounts for omitted variable
bias as in many applications, the implicit assumption that all confounding
variables are included is not easy to attain. Confidence intervals for
$GPG_{Qd}$ and $GPG_{Ud}$ are obtained for each area.

Furthermore, to validate the methodology proposal, a simulation study is
conducted. It has a particular design that mimics the generation of data of
the real application at hand. The estimators obtained with the novel approach
are compared with those using the OB approach, showing that the formers
outperform the latter.

We organize the remainder of the paper as follows. Section 2 revises the OB
approach and SAE models, while Section 3 details the novel proposal. The
estimation of the gender gap by economic activities in the region of Galicia
(Spain) is presented in Section 4 and the simulation study in Section 5.
Finally, the main conclusions and a brief discussion of future research is
given in Section 6.

##  2 Background

We use the notation which is usually considered for the OB decomposition and
SAE models. The beta parameter appears in both formulations, but with a
different role in each one.

###  2.1 Pay gap decomposition

Let $Y_{g}=log(W_{g})$ be the natural log of hourly earnings for $g=m,w$; men
and women, respectively. And let $\mathbf{X}_{g}$ ;$(X_{g1},\ldots,X_{gp})$
represent the set of available explanatory variables. Typically, these
explanatory variables are the education level, the age of the worker or the
work experience.

Assume that samples of sizes $n_{m},n_{w}$ are available for men and women,
respectively, and the following models:

$Y_{gi}=\mathbf{X}_{gi}\boldsymbol{\beta}_{g}+\epsilon_{gi},\quad
i=1,...,n_{g},\quad g\in(m,w);$

where $\boldsymbol{\beta}_{g}=(\beta_{g1},\ldots,\beta_{gp})^{\prime}$ is an
unknown vector of parameters and $\epsilon_{gi}$ is independent, with
$E(\epsilon_{g})=0$.

Then,

$\Delta=E({Y}_{m})-E({Y}_{w})=\overline{\mathbf{X}}_{m}\boldsymbol{\beta}_{m}-\overline{\mathbf{X}}_{w}\boldsymbol{\beta}_{w},$

where, $\bar{\mathbf{X}}_{m}$ and $\bar{\mathbf{X}}_{w}$ are the corresponding
sample means.

Now, let $\beta_{r}$ be the auxiliary coefficient; then, adding and
subtracting $\overline{\mathbf{X}}_{m}\boldsymbol{\beta}_{r}$ and
$\overline{\mathbf{X}}_{w}\boldsymbol{\beta}_{r}$, the above difference can be
written as:

$\Delta=(\overline{\mathbf{X}}_{m}-\overline{\mathbf{X}}_{w})\boldsymbol{\beta}_{r}+\overline{\mathbf{X}}_{m}(\boldsymbol{\beta}_{m}-\boldsymbol{\beta}_{r})+\overline{\mathbf{X}}_{w}(\boldsymbol{\beta}_{r}-\boldsymbol{\beta}_{w}).$

Frequently, the discrimination is directed toward one of the groups (women, in
our case), which determines the value of the auxiliary coefficient, so that
$\beta_{r}=\beta_{m}$ and then,

$\Delta=(\overline{\mathbf{X}}_{m}-\overline{\mathbf{X}}_{w})\boldsymbol{\beta}_{m}+\overline{\mathbf{X}}_{w}(\boldsymbol{\beta}_{m}-\boldsymbol{\beta}_{w}).$

where

$Q=(\overline{\mathbf{X}}_{m}-\overline{\mathbf{X}}_{w})\boldsymbol{\beta}_{m},$

is the component that is explained by group differences in the predictors (the
“quantify effect” or “the explained part”) and

$U=\overline{\mathbf{X}}_{w}(\boldsymbol{\beta}_{m}-\boldsymbol{\beta}_{w}).$

is the unexplained component, usually attributed to discrimination. Notice
that $U$ also captures all the potential effects of differences in unobserved
variables.

The estimation of the decomposition is straightforward. Let
$\hat{\boldsymbol{\beta}}_{m}$ and $\hat{\boldsymbol{\beta}}_{w}$ be the
least-squares estimates for $\boldsymbol{\beta}_{m}$ and
$\boldsymbol{\beta}_{w}$ obtained separately from the two samples. Then,

$\displaystyle\begin{array}[]{lll}\hat{\Delta}&=&{\bar{Y}}_{m}-{\bar{Y}}_{w}=\hat{Q}+\hat{U},\\\
\hat{Q}&=&(\bar{\mathbf{X}}_{m}-\bar{\mathbf{X}}_{w})\hat{\boldsymbol{\beta}}_{m},\\\
\hat{U}&=&\bar{\mathbf{X}}_{w}^{{}^{\prime}}(\hat{\boldsymbol{\beta}}_{m}-\hat{\boldsymbol{\beta}}_{w}).\end{array}$

More details and other decompositions can be seen in Jann (2008) and Hlavac
(2018).

###  2.2 SAE models

The wage estimation by small areas is a typical SAE problem as it is assumed
that the sample sizes are too low in many areas. The most extended approach in
SAE problems is to give a model-based estimator that usually includes random
effects due to the areas. To clarify the presentation we have only considered
linear models, as these are widely used in SAE applications and are well
suited to analyzing the problem at hand. Nevertheless, the proposal in the
paper can be adapted to consider other additive non-linear models as those
which include monotone and spline models such as those in Opsomer et al.
(2008) and Rueda and Lombardía (2012).

Specifically, the so-called nested error regression model with sampling
weights, which, for the $i$-th worker in area $d$, is defined as follows:

$Y_{di}=\mathbf{X}_{di}\boldsymbol{\beta}_{d}+\mathbf{Z}_{di}\mathbf{u}_{d}+\epsilon_{di},d=1,\ldots,D,i=1,\ldots,n_{d};$
(3)

where $Y_{di}$ is the (log) wage, $\mathbf{X}_{di}=(X_{1di},\ldots,X_{pdi})$
is a $p$ vector of explanatory variables, $\boldsymbol{\beta}_{d}$ is the
vector of unknown parameters; $\mathbf{Z}_{di}\mathbf{u}_{d}$ are the random
effects, with $\mathbf{Z}_{di}$ being a $q$ dimensional vector also defined
from the explanatory variables, and $\mathbf{u}_{d}\in
N(0,\boldsymbol{\Sigma}_{u})$, where $\boldsymbol{\Sigma}_{u}$ is a matrix of
dimension $q$. $\mathbf{u}_{d}$ is independent of the model error
$\epsilon_{di}$, which are assumed to be independent
$N(0,\sigma_{e}^{2}w_{di}^{-1})$, where $w_{di}$ is the sampling weight for
individual $i$ in area $d$. In these models the covariance of the response
variable in area $d$ is
$\mathbf{V}_{d}=\mathbf{Z}_{d}\boldsymbol{\Sigma}_{u}\mathbf{Z}_{d}+\boldsymbol{\Sigma}_{e}$
where $\boldsymbol{\Sigma}_{e}=\sigma_{e}^{2}{W_{d}}^{-1}$ and
${W_{d}}=diag(w_{di})_{n_{d}\times n_{d}}$. Two particular interesting cases
are considered by Battese et al. (1988) and Dempster et al. (1981).

In SAE applications, researchers are mainly interested in deriving estimators
for the response variable in small areas; they usually work with the
conditional means as the parameters of interest, which are defined as follows:

$\mu_{d}=E(Y_{d}|\mathbf{u}_{d})=\bar{\mathbf{X}}_{d}\boldsymbol{\beta}_{d}+\bar{\mathbf{Z}}_{d}\mathbf{u}_{d},d=1,\ldots,D$

The estimation process starts with the variance components
$(\boldsymbol{\Sigma}_{u},\boldsymbol{\Sigma}_{e})$, which are estimated by
the residual maximum likelihood (REML) method; the empirical best linear
unbiased estimators of $\boldsymbol{\beta}_{d}$ and the empirical best linear
unbiased predictors of $\mathbf{u}_{d}$ are then obtained. Furthermore,
estimators for $\mu_{d}$ are derived as follows:

$\hat{\mu}_{d}=\bar{\mathbf{X}}_{d}\hat{\boldsymbol{\beta}}_{d}+\bar{\mathbf{Z}}_{d}\hat{\mathbf{u}}_{d},d=1,\ldots,D$

For details on the estimation process and further learning about SAE methods,
the reader is referred to the monographs of Jiang (2007) and Rao and Molina
(2015).

##  3 The novel proposal

Assume that we have wage data for men and women and explanatory variables from
samples from $D$ small areas, and that we are interested in deriving estimates
for the two components of the $GPG$.

A three-step approach is proposed to derive GPG estimators. The first step is
a model selection step; in the second step, preliminary area-specific
estimators for $GPG_{Q}$ and $GPG_{U}$ are obtained from a GPG decomposition
adapted to small areas; and finally, in the third step, a Monte-carlo
algorithm is consider to derive corrected bias $GPG_{Q}$ and $GPG_{U}$
estimators and confidence intervals.

First Step: Model Selection

There are different approaches to model selection, and the task is complicated
when random effects are present. The discussions by Lombardia et al. (2017),
where the GAIC is first presented, and by Fan and Li (2012), where iterative
regularization methods are proposed, are very instructive. Besides proposing a
method for making the model selection, a list of candidate models must be
defined. From the methodological point of view, it can be as large as the
practitioner decides. From the computational and interpretative points of
view, a reduced list of candidate models are recommended.

The strategy adopted in this paper is to define a reduced list of candidate
models and then to use the GAIC statistic to select a model from this list. In
particular, for applications with a reduced number and well-known, explanatory
variables, we recommend the use of prior knowledge to decide an initial
selection of fixed and random effects. Typically, in small area applications,
the variable determining the areas is modelled using random effects. Other
candidate models may be defined by considering the interactions of that
variable, with other explanatory variables.

The group (men/women) with the larger sample size is used to select the model.
In $GPG$ applications, it is usually that of men ($m$). The models from which
we make the choice differ in the explanatory variables selected and the type
of effects (fixed or random) that describe the model and are formulated as the
general model (3), as follows:

$M:Y_{mdi}=\mathbf{X}_{mdi}\boldsymbol{\beta}_{md}+\mathbf{Z}_{mdi}\mathbf{u}_{md}+\epsilon_{mdi},\quad\
d=1,\ldots,D,i=1,\ldots,n_{d}.$

The selection of the best model is done using an $AIC$ statistic (Akaike,
1973); specifically, we use an adapted version for unit level models of the
$xGAIC$ proposed by Lombardía et al. (2017) with the following expression for
the model $M$

$xGAIC=-2log(l_{x}(M))+xGDF$

where $l_{x}(M)$ is the quasi-loglikelihood of the model $M$, which considers
the focus on the random effect and the total variability and $xGDF$ is a
measure of the complexity of the model, as follows:

$\log(l_{x}(M))=-\frac{1}{2}D\log(2\pi)-\frac{1}{2}\log|\mathbf{V}_{Y_{m}}|-\frac{1}{2}(\mathbf{Y}_{m}-\boldsymbol{\mu}_{m})^{{}^{\prime}}\mathbf{V}_{Y_{m}}^{-1}(\mathbf{Y}_{m}-\boldsymbol{\mu}_{m})$

and

$xGDF=\sum_{d=1}^{D}\frac{\partial
E(\widehat{\boldsymbol{\mu}}_{md})}{\partial\boldsymbol{\mu}_{md}}=\sum_{d=1}^{D}\sum_{i=1}^{n_{md}}\sum_{j=1}^{n_{md}}V_{mdij}^{-1}Cov(\widehat{\mu}_{mdi},Y_{mdj})$

where $\mathbf{Y}_{m}$ is the $\sum_{d=1}^{D}n_{d}\times 1$ vector with
elements $Y_{mdi}$, $\boldsymbol{\mu}_{m}$ is the vector of the conditional
means $\boldsymbol{\mu}_{m}=E(\mathbf{Y}_{m}|\mathbf{u})$ with elements
$\boldsymbol{\mu}_{mdi}=\mathbf{X}_{mdi}\boldsymbol{\beta}_{md}+\mathbf{Z}_{mdi}\mathbf{u}_{md}$,
$\mathbf{V}_{Y_{m}}$ is the conditional covariance matrix, which is diagonal
with elements $\mathbf{V}_{md}$ according to Section 2.2 and $V_{mdij}^{-1}$
is the $ij$-element of $\mathbf{V}_{md}^{-1}$.

Parametric bootstrap is used to estimate $xGDF$, as the analytic values are
difficult to obtain in a similar way to Lombardía et al. (2017). The model
with the lowest value of $xGAIC$ is selected.

Second Step: GPG decomposition. Preliminary estimators.

Let $X_{m},Z_{m}$ and $X_{w}$,$Z_{w}$ be the design matrices representing the
explanatory variables and the random effects on the selected model for men and
women, respectively. For a given area, $d$, consider
$\Delta_{d}=\mu_{md}-\mu_{wd}$, which is decomposed as follows:

$\Delta_{d}=\mu_{md}-\mu_{wd}=\bar{\mathbf{X}}_{md}\boldsymbol{\beta}_{md}+\bar{\mathbf{Z}}_{md}\mathbf{u}_{md}-\bar{\mathbf{X}}_{wd}\boldsymbol{\beta}_{wd}-\bar{\mathbf{Z}}_{wd}\mathbf{u}_{wd}.$

$\displaystyle\begin{array}[]{lll}\Delta_{d}&=&Q_{d}+U_{d}\\\
{Q}_{d}&=&(\bar{\mathbf{X}}_{md}-\bar{\mathbf{X}}_{wd})\boldsymbol{\beta}_{md}+(\bar{\mathbf{Z}}_{md}-\bar{\mathbf{Z}}_{wd})\mathbf{u}_{md},\\\
{U}_{d}&=&\bar{\mathbf{X}}_{wd}(\boldsymbol{\beta}_{md}-\boldsymbol{\beta}_{wd})+\bar{\mathbf{Z}}_{wd}(\mathbf{u}_{md}-\mathbf{u}_{wd}).\end{array}$

and,

${GPG}_{Qd}={GPG}_{d}\frac{Q_{d}}{\Delta_{d}}\ \quad\ \mbox{and}\quad\
{GPG}_{Ud}={GPG}_{d}\frac{{U}_{d}}{\Delta_{d}}.$

where

$GPG_{d}=\frac{E(W_{md})-E(W_{wd})}{E(W_{md})},$

Small area estimators for $Y_{gd}$, $\Delta_{d}$, $U_{d}$ and $Q_{d}$ (as in
Rao and Molina, 2015) and estimators for $E(W_{gd}),g\in(m,w)$ $d=1,...,D$,
are defined as follows:

$\displaystyle\begin{array}[]{lll}\hat{\Delta}_{d}&=&\bar{\mathbf{X}}_{md}\hat{\boldsymbol{\beta}}_{md}+\bar{\mathbf{Z}}_{md}\hat{\mathbf{u}}_{md}-\bar{\mathbf{X}}_{wd}\hat{\boldsymbol{\beta}}_{wd}-\bar{\mathbf{Z}}_{wd}\hat{\mathbf{u}}_{wd},\\\
\hat{Q}_{d}&=&(\bar{\mathbf{X}}_{md}-\bar{\mathbf{X}}_{wd})\hat{\boldsymbol{\beta}}_{md}+(\bar{\mathbf{Z}}_{md}-\bar{\mathbf{Z}}_{wd})\hat{\mathbf{u}}_{md},\\\
\hat{U}_{d}&=&\hat{\Delta}_{d}-\hat{Q}_{d}.\end{array}$ (8)

$\widehat{E(W_{gd})}=\overline{exp(\hat{Y}_{gd})},g\in(m,w),d=1,...,D;$ (9)

From (9), we have that:

$\widehat{GPG_{d}}=\frac{\widehat{E(W_{md})}-\widehat{E(W_{wd})}}{\widehat{E(W_{md})}},$
(10)

Third Step: bias corrected estimators and confidence intervals.

Omitted variable bias is taken into account as in many applications, the
implicit assumption that all confounding variables are included is not easy to
attain. The bias in $Q_{d}$ is derived by splitting the men sample into two
parts I times. The value of $B_{d}^{i}$, for iteration $i=1,...,I$, is
expected to be zero when comparing the populations representing these two
parts. A non-zero estimated value is assumed to be due to the omitted
variables and is denoted by $\hat{B}^{i}_{d}$, for iteration $i$. From these
values, the bias in $Q_{d}$ is estimated, and bias-corrected estimators for
$GPG_{Qd}$ and $GPG_{Ud}$ are derived, as explained below. Iterations with no
men or women in any of the activities are discarded.

For each iteration $i=1,\ldots,I$ repeat:

  1. 1.

Randomly split the data $S$ in two parts of size $\frac{n}{2}$: $S_{1}^{i}$
and $S_{2}^{i}$. The individuals of $S_{1}^{i}$ are selected randomly, in such
a way that those that have not been selected form part of $S_{2}^{i}$.
Therefore, the union of these two sets is the complete sample. In the
following, subindex 1 or 2 indicates the data coming from $S_{1}^{i}$ or
$S_{2}^{i}$, respectively.

  2. 2.

The model selected is fitted to men and women’s data in $S_{1}^{i}$, and
$\hat{Q}_{1d}^{i}$, $\hat{\Delta}_{1d}^{i}$, and $\widehat{GPG}_{1d}^{i}$ are
calculated according to (8) and (10).

  3. 3.

For each area $d$ and iteration $i$, the bias is calculated considering men’s
data from $S_{1}^{i}$ and $S_{2}^{i}$, as follows:

$\hat{B}_{d}^{i}=(\bar{\mathbf{X}}_{1md}-\bar{\mathbf{X}}_{2md})\hat{\boldsymbol{\beta}}_{1md}+(\bar{\mathbf{Z}}_{1md}-\bar{\mathbf{Z}}_{2md})\hat{u}_{1md}.$

Then, the final bias term is estimated as:

$\hat{B}_{d}=\frac{1}{I}\sum_{i=1}^{I}\hat{B}_{d}^{i},$

and the corrected bias estimates of ${Q}_{d}$ and ${U}_{d}$ are:

$\hat{Q}_{d}^{B}=\hat{Q}_{d}+\hat{B}_{d},$
$\hat{U}_{d}^{B}=\hat{\Delta}_{d}-\hat{Q}_{d}^{B},$

where $\hat{Q}_{d}$ and $\hat{\Delta}_{d}$ are from the complete sample, as in
(8).

Now, the corrected bias estimates of the explained and unexplained parts of
$GPG$ are derived from the latter and from (10) as:

$\widehat{GPG}_{Qd}=\widehat{GPG}_{d}\frac{\hat{Q}_{d}^{B}}{{\hat{\Delta}_{d}}}\
\quad\ \mbox{and}\quad\
\widehat{GPG}_{Ud}=\widehat{GPG}_{d}\frac{\hat{U}_{d}^{B}}{{\hat{\Delta}_{d}}}$
(11)

Moreover, confidence intervals are also derived as follows:

From $S_{1}^{i}$, $i=1,\ldots,I$, and (9), we have that:

$\hat{Q}_{d}^{Bi}=\hat{Q}_{1d}^{i}+\hat{B}_{d},\quad\mbox{and}\quad\hat{U}_{d}^{Bi}=\hat{\Delta}_{1d}^{i}-\hat{Q}_{d}^{Bi}.$

$\widehat{GPG}_{1d}^{i}=\frac{\widehat{E(W_{1md})}-\widehat{E(W_{1wd})}}{\widehat{E(W_{1md})}},$
$\widehat{GPG}_{Qd}^{i}=\widehat{GPG}_{1d}^{i}\frac{\hat{Q}_{d}^{Bi}}{\hat{\Delta}_{1d}^{i}},$
$\widehat{GPG}_{Ud}^{i}=\widehat{GPG}_{1d}^{i}\frac{\hat{U}_{d}^{Bi}}{\hat{\Delta}_{1d}^{i}},$

Now, let us define:

$\displaystyle\begin{array}[]{lll}V(\widehat{GPG}_{Qd})&=&\frac{1}{I}\sum_{i=1}^{I}\left(\widehat{GPG}_{Qd}^{i}-\widehat{\overline{GPG}}_{Qd}^{B}\right)^{2},\\\
V(\widehat{GPG}_{Ud})&=&\frac{1}{I}\sum_{i=1}^{I}\left(\widehat{GPG}_{Ud}^{i}-\widehat{\overline{GPG}}_{Ud}^{B}\right)^{2},\end{array}$

and,

$\widehat{\overline{GPG}}_{Qd}^{B}=\frac{1}{I}\sum_{i=1}^{I}\widehat{GPG}_{Qd}^{i}.$
$\widehat{\overline{GPG}}_{Ud}^{B}=\frac{1}{I}\sum_{i=1}^{I}\widehat{GPG}_{Ud}^{i}.$

while the confidence intervals are derived as follows:

$\displaystyle{}\begin{array}[]{lll}\widehat{{GPG}}_{Qd}\pm
Z_{\alpha/2}\sqrt{V(\widehat{GPG}_{Qd})}\\\ \widehat{{GPG}}_{Ud}\pm
Z_{\alpha/2}\sqrt{V(\widehat{GPG}_{Ud})}\end{array}$ (15)

The aim of this document is to provide an easy-to-apply methodology for
constructing confidence intervals to allow for inferences to be made. The
post-selection inference is an interesting issue that could be considered in
the proposed methodology following the ideas of Berk et al. (2013), Lee et al.
(2016) and Tibshirani et al. (2016), among others; but it is outside the scope
of this document.

##  4 Estimation of the GPG by ecomomic activities in Galicia.

The data come from the Structure of Earnings Survey (SES) conducted by the
Spanish National Institute of Statistics (INE) in 2014. The subset considering
in this paper, corresponds to the Galicia region at the Northwest of Spain.
The survey is conducted in a similar way to wage structure surveys in other
European countries; it uses a two-level sampling of national companies, a
stratified sampling in the first stage for local units and systematic sampling
to select workers at those units. The economic activities are defined by the
Statistical Classification of Economic Activities in the European Community
(NACE09 Rev.2) at the division level.

Galicia has an economy strongly linked to natural resources, and its labor
market has experienced growth in recent years. The $GPG$, it has been reduced
from $15.5\%$ in 2010 to $13.5\%$ in 2016. This latter value is below that of
the EU ($16.3\%$), and below that of Spain ($15.1\%$) in that year.

The survey collects information on non-self-employed workers who work in
establishments with at least 10 employees and covers a wide range of private
sectors (industry, construction, commerce, catering business, transport,
financial intermediation, …) excluding the primary sector. The total sample
size in Galicia being 10276. The list of explanatory variables considered in
this study are included in Table 1, the selection being similar to that in
Aláez et al., (2000, 2001, 2003), De la Rica et al. (2008) or Moral-Arce et
al. (2011).

The economic activities, which are defined by the Statistical Classification
of Economic Activities in the European Community (NACE09 Rev.2) at the
division level, are the areas of interest. Tables 10 and 11 in the Appendix
include the list of codes and the description of the NACE09 at the division
level, excluding the primary sector. In the present study, we have information
for $D=78$ of the 84 NACE09’s. SAE methods are appealing because the sample
sizes in some activities are very small. Specifically, the quartiles of the
sample sizes distribution in the areas are $Q_{1}=26$, $Q_{2}=50$ and
$Q_{3}=104$, for men, and $Q_{1}=15$, $Q_{2}=26$ and $Q_{3}=50$, for women.

Table 2 gives the main statistics for men and women regarding labor
characteristics. The numbers in the table show a greater percentage of women
in higher education levels, 41% against 33% of men, and that almost half the
men (48%) are employed in the manufacturing industry or in public
administration, while public administration (41%) and the retail trade (16%)
are the activities with higher percentages for women. Regarding their type of
work, women are more frequently on part-time work (27% vs 10% for men).

Name | Variable | Description
---|---|---
Wage per hour | Y | Gross hourly earnings from employment
Gender |  | Men, Women
Experience | $X_{1}$ | Number of years in the actual enterprise
Age | $X_{2}$ | Years of the employee
Education |  | (Primary)
| $X_{3}$ | Secondary
| $X_{4}$ | Higher
Occupation |  | (Professionals and managers)
| $X_{5}$ | Tecnic
| $X_{6}$ | Operators
| $X_{7}$ | Services and sales workers
| $X_{8}$ | Non skilled workers
Contract | $X_{9}$ | (Long term), Short term
Type of work | $X_{10}$ | (Full time), Part time
Size of the enterprise |  | (Between 10 to 19 employees)
| $X_{11}$ | Between 20 to 49 employees
| $X_{12}$ | Between 50 to 99 employees
| $X_{13}$ | Between 100 to 199 employees
| $X_{14}$ | Between 200 to 499 employees
| $X_{15}$ | More than 500 workers
Market |  | (Local or regional)
| $X_{16}$ | National
| $X_{17}$ | UE
| $X_{18}$ | International
Aggregated economic activity |  | (Energy)
| $X_{19}$ | Manufacturing industry
| $X_{20}$ | Construction
| $X_{21}$ | Retail trade
| $X_{22}$ | Transportation, storage and Accommodation
| $X_{23}$ | Information and Communication
| $X_{24}$ | Finance and insurance
| $X_{25}$ | Professionals
| $X_{26}$ | Public Administration
| $X_{27}$ | Other service
Economic activity for two digits | $X_{28}-X_{104}$ | (7) NACE09 Rev. 2, division level, two digits
Table 1:  Variables in the survey, in brackets the category used as a reference. For the variable Economic activity for two digits, the reference category is the activity 7, because in Galicia there are not people working in the activity 6 (see Appendix).  | Men | Women
---|---|---
Experience (mean years) |  |
| 11.16 | 10.43
Education (%) |  |
Primary | 0.19 | 0.13
Secondary | 0.48 | 0.46
Higher | 0.33 | 0.41
Contract (%) |  |
Long-term | 0.77 | 0.77
Short-term | 0.23 | 0.23
Type of work (%) |  |
Full time | 0.90 | 0.73
Part time | 0.10 | 0.27
Economic activity (%) |  |
Energy | 0.03 | 0.01
Manufacturing industry | 0.27 | 0.13
Construction | 0.09 | 0.01
Retail trade | 0.13 | 0.16
Transportation and storage and Accommodation | 0.08 | 0.07
Information and Communication | 0.03 | 0.02
Finance and insurance | 0.03 | 0.03
Professionals | 0.10 | 0.14
Public Administration | 0.21 | 0.41
Other services | 0.02 | 0.02
Table 2:  Summary of labor force characteristics of the people in the survey
by gender.

The distribution of the response variable ($Y$), the wage per hour in
logarithm, is shown in Figure 1 for men and women; the shape of both curves
being quite similar, though that of men is shifted to the right, which
indicates an increase in salaries. Specifically, the quartiles are Q1 = 7.6,
Q2 = 9.6 and Q3 = 13.6 for men and Q1 = 6.4, Q2 = 8.1 and Q3 = 12.1 for women.
Finally, Figure 2 shows important differences in mean wages per hour across
activities and between sexes. These differences range from 9 euros in the
median of the activity retail trade and the maximum of 18 euros for those
employed in finance and insurance activities.

Figure 1:  Density of ln(wage/hour) by sex in Galicia. Figure 2:  Wage per
hour by sex and economic activity in Galicia. The codes for the activities
are: 2 Energy, 3 Manufacturing industry, 4 Construction, 5 Retail trade, 6
Transportation and storage and Accommodation, 7 Information and Communication,
8 Finance and insurance, 9 Professionals, 10 Public Administration and 11
Other services.

The candidate list of models have been built as follows; First, the
explanatory variables more often used in similar studies have been selected
and have been described in Table 1. Next, the candidate models are defined
combining fixed effect for Experience, Age, Education, Occupation, Contract,
Type of work, Size of the enterprise, and Market with the effects by the
Economic Activities. The latter can be incorporated in the model using fixed
effects for aggregated levels ($X19$ to $X27$), or fixed effects for each
activity ($X28$ to $X104$), or a random effect. Moreover, when a random effect
is used to model the activities, the interactions with Experience, Education,
and Occupation may also be considered.

Regarding the notation, $u_{d}$ is the random effect for the Economic
Activities, $u_{dj},j=1,2$ for the interaction with Education, $u_{dk}$
$k=1,2,3,4$ for the interaction with Occupation; and $v_{d}$ for the
interaction with Experience. Only a selection of the most interesting models
is presented in Table 3 to simplify the exposition. Specifically, eight models
labeled as $M1$ to $M8$ are presented, models $M6$ and $M7$ are defined with
fixed effects, and the rest with a combination of fixed and random effects.
Models defined by other combinations of effects have been considered but have
been discarded because the results are far from being better than those
obtained with other models.

Label | Fixed effects | Random Effects | Fixed parameters (nº) | Random parameters (nº)
---|---|---|---|---
$M1$ |  $X_{1}$-$X_{18}$ | $u_{d}$ | 18 | 1
$M2$ |  $X_{1}$-$X_{18}$ | $v_{d}$ | 18 | 1
$M3$ |  $X_{1}$-$X_{18}$ |  $u_{d}$, $v_{d}$ | 18 | 2
$M4$ |  $X_{1}$-$X_{18}$ | $u_{dj}$ | 18 | 2
$M5$ |  $X_{1}$-$X_{18}$ | $u_{dk}$ | 18 | 4
$M6$ |  $X_{1}$-$X_{27}$ | - | 27 | 0
$M7$ |  $X_{1}$-$X_{18}$, $X_{28}$-$X_{104}$ | - | 95 | 0
$M8$ |  $X_{1}$-$X_{18}$, $X_{28}$-$X_{104}$ | $v_{d}$ | 95 | 1
Table 3:  Models used to fit the data.

In order to explore differences between models, the estimated bias for each
activity has been calculated for each candidate model and is shown in Figure
3. The activities are sorted in increasing order of bias according to the $M1$
model. Figure 3 shows important bias in some activities: $88=$ Social work
activities without accommodation, $90=$ Creative, arts and entertainment
activities or $39=$ Remediation activities and other waste management
services, among others. Moreover, the bias corrected estimator of the
explained and unexplained parts of the $GPG$ for each economic activity and
model, as shown in the expression (11), are shown in Figures 4 and 5,
respectively, sorting the activities by the values of $\widehat{GPG}_{Qd}$ and
$\widehat{GPG}_{Ud}$ for $M1$. These Figures show that while estimates for
$\widehat{GPG}_{Qd}$ are similar between models, there are important
differences in $\widehat{GPG}_{Ud}$ across models.

Figure 3: Bias for $Q_{d}$ ($\hat{B}_{d}$), for each activity $d$ defined by
the NACE09. Figure 4:  $\widehat{GPG}_{Qd}$ for all the models considered and
for each activity $d$ defined by the NACE09. Figure 5:  $\widehat{GPG}_{Ud}$
for all the models considered and for each activity $d$ defined by the NACE09.

The model selected, using the proposal described in Section 3, is $M1$. The
parameter estimates for $M1$ are shown in Table 4, which explains that the
wage of a worker increase with higher education, working in an enterprise with
a national, UE or international market, and belonging to a large enterprise.
The coefficients related to occupations show negative signs as Professionals
and managers is the reference category. The results are consistent with those
from other studies (Aláez et al. (2003), De la Rica (2008) and Moral-Arce et
al. (2011)).

| Estimate | Std. Error | p-value
---|---|---|---
Fixed effects |  |  |
(Intercept) | 2.299 | 0.033 | 0.000
Experience ($X_{1}$) | 0.007 | 0.001 | 0.000
Age ($X_{2}$) | 0.005 | 0.000 | 0.000
Secondary ($X_{3}$) | 0.015 | 0.011 | 0.188
Higher ($X_{4}$) | 0.122 | 0.014 | 0.000
Technic ($X_{5}$) | -0.248 | 0.015 | 0.000
Operators ($X_{6}$) | -0.439 | 0.016 | 0.000
Services and sales workers ($X_{7}$) | -0.504 | 0.019 | 0.000
Non-skilled workers ($X_{8}$) | -0.527 | 0.016 | 0.000
Short term contract ($X_{9}$) | -0.098 | 0.010 | 0.000
Part time work ($X_{10}$) | -0.017 | 0.013 | 0.210
Between 20 to 49 workers ($X_{11}$) | -0.029 | 0.012 | 0.028
Between 50 to 99 workers ($X_{12}$) | 0.047 | 0.014 | 0.003
Between 100 to 199 workers ($X_{13}$) | 0.042 | 0.014 | 0.011
Between 200 to 499 workers ($X_{14}$) | 0.081 | 0.014 | 0.000
More than 500 workers ($X_{15}$) | 0.113 | 0.016 | 0.000
National Market ($X_{17}$) | 0.052 | 0.010 | 0.000
UE Market ($X_{18}$) | 0.076 | 0.018 | 0.001
International Market ($X_{19}$) | 0.072 | 0.016 | 0.000
Random effects |  |  |
$\sigma_{u}^{2}$ | 0.021 |  |
Table 4:  Coefficient estimates for men and for $M3$.

Pay gap estimators have been obtained using the model $M1$ and the proposal in
Section 3. In Galicia, the pay gap between men and women ($\widehat{GPG}$) is
estimated to be of 11.3%, where 2.8% is due to the observed quantify effect
($\widehat{GPG_{Q}}$) and 8.5% to discrimination ($\widehat{GPG_{U}}$), as
Table 5 shows.

| Estimation | CI 95%
---|---|---
$\widehat{GPG_{Q}}$ | 0.028 | (-0.018,0.07)
$\widehat{GPG_{U}}$ | 0.085 | (0.057,0.14)
$\widehat{GPG}$ | 0.113 | -
Table 5:  Decomposition of gender pay gap in Galicia.

Moreover, Figures 6 and 7 show $\widehat{GPG}_{Qd}$ and $\widehat{GPG}_{Ud}$,
respectively, for each activity, $d$, with their confidence intervals.
Activities are ordered by the increasing value of the estimator.

It is worth commenting on the results from activities 84 and 85, which
correspond to the Public Administration and Education. They are useful for
validating the procedure because no discrimination is expected in these cases
as the Administration control the wages (Aláez et al., 2000, Moral-Arce et
al., 2011). Figure 6 shows no significant differences between the
characteristics of men and women, for activities 84 and 85; while Figure 7
shows, for 84 and 85, no significative differences from those of Galicia. The
results are consistent with what would be expected.

On the other hand, the results show very high discrimination in eight
activities: 30, 28, 19, 17, 42 , 33, 29 and 10, almost all of them are from
the industrial sector (see Figure 7); while differences between the
characteristics due to gender are not significantly higher for these
activities than for Galicia (Figure 6).

Figure 6:  $\widehat{GPG}_{Qd}$ with their confidence interval (95%) and for
each activity $d$ defined by the NACE09. The horizontal continuous line is the
$\widehat{GPG}_{Q}$ for Galicia. Figure 7:  Discrimination
($\widehat{GPG}_{Ud}$) with their confidence interval (95%) and for each
activity $d$ defined by the NACE09. The horizontal continuous line is the
$\widehat{GPG}_{U}$ for Galicia.

Taking as reference the $GPG$ of the region of Galicia, Table 9 shows the
activities with higher $\widehat{GPG_{d}}$ than Galicia, sorted from highest
to lowest. The activity at the top is Electricity, gas, steam and air
conditioning supply and it is also one of the activities where there are more
discrimination. The second activity at the top is Manufacture of wearing
apparel, and this is an example that high wage differences are not necessarily
accompanied by discrimination. This activity is very important in the Galician
business sector and employs approximately 8,600 workers, $1\%$ of the total
number of employed people in Galicia. On the other hand, note that only ten
activities out of 74 have positive discrimination significantly higher than
the region’s mean. These activities employ $8.2\%$ of the workers in the
region. It is also interesting to note that there are two activities in which
discrimination is significantly negative (in favor of women). These are
Warehousing and support activities for transportation and  Security and
investigation activities that belong to the service sector where women have a
more significant presence than men.

##  5 Simulations

The aim of this section is to validate the novel proposal and compare it with
the classical approach of OB. Moreover, the estimation error that results when
important explanatory variables are omitted, is also illustrated.

The simulated experiment imitates SES real data. The 30 activities which have
the largest sample sizes and the explanatory variables in Table 1 are
considered. The sample sizes take values from 21 to 111 for men and 7 to 82
for women. The model used to generate the data is defined as follows:

$\displaystyle\begin{array}[]{lll}Y_{gdi}&=&X_{1gdi}\beta_{1gd}+X_{3gd}\beta_{2gdi}+X_{4dg}\beta_{3gd}+X_{5gdi}\beta_{4gd}+X_{6gdi}\beta_{5gd}+X_{7gdi}\beta_{6gd}+X_{8gdi}\beta_{7gd}+\epsilon_{gdi}\\\
&&g\in\\{m,w\\};\quad\ d=1,...,30;\quad\ i=1,\ldots,n_{gd};\end{array}$ (18)

where $\epsilon_{gdi}\sim N(0,\sigma_{gd}^{2})$, with $\sigma_{gd}^{2}=0.1$
and $\beta_{jgd}$, $j=1,\ldots,7$ are the estimated values using model $M6$ .

Note that model (18) is much more complex than any of the models in the
candidate list. The generating model includes a fixed parameter for each
activity and explicative variable, giving a model with 8*30=240 parameters;
being the candidate models much more simple.

For each $d$, $GPG_{d}$, $GPG_{Qd}$ and $GPG_{Ud}$ are obtained using (1), (2)
and (5) and are assumed to be true values.

The models in Table 6 combine fixed and random effects, as in the real case.
Estimators for GPG components and each model are obtained using proposal in
Section 3.

Label | Fixed effects | Random Effects | Fixed parameters (nº) | Random parameters (nº)
---|---|---|---|---
$MS1$ |  $X_{1}$,$X_{3}$-$X_{8}$ | $u_{d}$ | 7 | 1
$MS2$ |  $X_{1}$,$X_{3}$-$X_{8}$ | $v_{d}$ | 7 | 1
$MS3$ |  $X_{1}$,$X_{3}$-$X_{8}$ |  $u_{d}$, $v_{d}$ | 7 | 2
$MS4$ |  $X_{1}$,$X_{3}$-$X_{8}$ | $u_{dj}$ | 7 | 2
$MS5$ |  $X_{1}$,$X_{3}$-$X_{8}$ | $u_{dk}$ | 7 | 4
$MS6$ |  $X_{1}$,$X_{3}$-$X_{8}$, $X_{19}$-$X_{27}$ | - | 16 | 0
$MS7$ |  $X_{1}$,$X_{3}$-$X_{8}$, $X_{28}$-$X_{56}$ | - | 36 | 0
$MS8$ |  $X_{1}$,$X_{3}$-$X_{8}$, $X_{28}$-$X_{56}$ | $v_{d}$ | 36 | 1
Table 6:  Models used to fit the data.

Now, the empirical mean square error (EMSE) of $\widehat{GPG}_{Qd}$ and
$\widehat{GPG}_{Ud}$ for all the models and for each activity $d$ are as
follows:

$EMSE(\widehat{GPG}_{Qd})=\frac{1}{100}\sum_{i=1}^{100}\left(\widehat{GPG}_{Qd}^{(i)}-{GPG}_{Qd}\right)^{2}$
$EMSE(\widehat{GPG}_{Ud})=\frac{1}{100}\sum_{i=1}^{100}\left(\widehat{GPG}_{Ud}^{(i)}-{GPG}_{Ud}\right)^{2}$

where $\widehat{GPG}_{Qd}^{(i)}$ and $\widehat{GPG}_{Ud}^{(i)}$ are the
estimations for each iteration $i=1,\ldots,100$.

Finally, we calculate the confidence intervals for $GPG_{Qd}$ and ${GPG_{Ud}}$
from (15) and give the percentage of times that the true values $GPG_{Qd}$ and
${GPG_{Ud}}$ are in the corresponding confidence interval.

Tables 7 and 8 present the results of the simulations, EMSE values and
coverage rates of confidence intervals, respectively. Table 7 shows the mean
of the EMSE of $\widehat{GPG}_{Qd}$ and $\widehat{GPG}_{Ud}$ obtained from the
models $MS1$ to $MS8$; the OB model and the model selected by xGAIC (XG) are
also included. The lefthand part of the table corresponds to the model with
all the explanatory variables and the righthand side shows the results when
the information relative to education ($X_{3}$ and $X_{4}$) is removed.

The best models, in terms of EMSE, are $MS5$ and $XG$. When relevant
information is missed out the error increases in all the models, but the
increase is higher for OB and the MS6 model. The conclusions are similar
regarding the coverage rates of 95% confidence intervals in Table 8. From
these results, it can be said that the proposed methodology is robust to the
potential effects of differences in unobserved variables; to be precise, much
more than the OB approach.

In addition, in Tables 7 and 8, you can see that the models selected by the
$xGAIC$ have a good behavior, because they are among those that have a lower
EMSE and a greater coverage. When the model is fitted with all the variables,
the $xGAIC$ selects the MS5 model $82\%$ of times, $15\%$ of times $MS4$ and
$3\%$ $MS8$. When $X_{3}$ and $X_{4}$ are simultaneously removed from the
working model, the $xGAIC$ selects the model $MS5$ $98\%$ of times and $2\%$
the model $MS8$.
Finally, this study highlights the gain of our methodology with respect to the
OB methodology, especially for estimating ${GPG}_{Ud}$.

|  |  | Without $X_{3}$, $X_{4}$ |
---|---|---|---|---
| $GPG_{Q}$ | $GPG_{U}$ | $GPG_{Q}$ | $GPG_{U}$
OB | 0.0049 | 0.0067 | 0.0060 | 0.0089
MS1 | 0.0046 | 0.0045 | 0.0054 | 0.0057
MS2 | 0.0057 | 0.0066 | 0.0067 | 0.0076
MS3 | 0.0060 | 0.0048 | 0.0066 | 0.0058
MS4 | 0.0045 | 0.0044 | 0.0054 | 0.0057
MS5 | 0.0044 | 0.0044 | 0.0052 | 0.0054
MS6 | 0.0049 | 0.0067 | 0.0061 | 0.0079
MS7 | 0.0046 | 0.0048 | 0.0053 | 0.0060
MS8 | 0.0060 | 0.0049 | 0.0065 | 0.0060
XG | 0.0045 | 0.0043 | 0.0052 | 0.0054
Table 7:  EMSE of $GPG_{Q}$ and $GPG_{U}$ for the model with and without $X_{3}$ and $X_{4}$. OB is the Oaxaca-Blinder model and XG is the model selected by the $xGAIC$ |  |  | Without $X_{3}$, $X_{4}$ |
---|---|---|---|---
| $GPG_{Q}$ | $GPG_{U}$ | $GPG_{Q}$ | $GPG_{U}$
OB | 86.4 | 55.0 | 80.1 | 43.7
MS1 | 86.4 | 78.9 | 80.5 | 73.7
MS2 | 83.6 | 57.4 | 79.0 | 54.6
MS3 | 88.7 | 81.4 | 84.2 | 77.8
MS4 | 90.1 | 83.0 | 80.5 | 73.7
MS5 | 90.8 | 84.5 | 85.1 | 81.5
MS6 | 86.4 | 55.0 | 80.2 | 45.9
MS7 | 86.4 | 81.7 | 80.2 | 76.2
MS8 | 88.7 | 83.2 | 83.9 | 79.4
XG | 90.7 | 84.4 | 85.0 | 81.0
Table 8:  Coverage of the $GPG_{Q}$ and $GPG_{U}$ confidence intervals for the model with and without $X_{3}$ and $X_{4}$. OX is the Oaxaca-Blinder model and XG is the model selected by the $xGAIC$ NACE09 | $\widehat{GPG}_{d}$ | $\widehat{GPG}_{Qd}$ | $\widehat{GPG}_{Qd}^{-}$ | $\widehat{GPG}_{Qd}^{+}$ | $\widehat{GPG}_{Ud}$ | $\widehat{GPG}_{Ud}^{-}$ | $\widehat{GPG}_{Ud}^{+}$
---|---|---|---|---|---|---|---
35 | 0.40 | 0.08 | -0.02 | 0.18 | 0.32 | 0.17 | 0.47
14 | 0.35 | 0.21 | 0.12 | 0.30 | 0.17 | -0.03 | 0.37
96 | 0.32 | 0.14 | 0.01 | 0.26 | 0.19 | 0.05 | 0.33
19 | 0.28 | -0.00 | -0.32 | 0.31 | 0.28 | 0.12 | 0.45
17 | 0.27 | 0.01 | -0.11 | 0.13 | 0.27 | 0.16 | 0.39
61 | 0.27 | 0.13 | 0.05 | 0.21 | 0.15 | 0.05 | 0.24
28 | 0.25 | -0.05 | -0.23 | 0.14 | 0.29 | 0.13 | 0.44
65 | 0.24 | 0.12 | -0.02 | 0.26 | 0.11 | -0.05 | 0.28
94 | 0.22 | 0.03 | -0.12 | 0.17 | 0.21 | -0.07 | 0.48
50 | 0.21 | 0.08 | -0.20 | 0.36 | 0.14 | -0.05 | 0.32
29 | 0.21 | 0.00 | -0.08 | 0.08 | 0.21 | 0.14 | 0.28
45 | 0.20 | 0.12 | -0.01 | 0.24 | 0.08 | -0.07 | 0.23
10 | 0.20 | -0.01 | -0.07 | 0.05 | 0.21 | 0.11 | 0.31
13 | 0.19 | 0.05 | -0.04 | 0.14 | 0.14 | 0.00 | 0.28
18 | 0.19 | 0.05 | -0.03 | 0.12 | 0.15 | 0.05 | 0.24
69 | 0.19 | 0.13 | -0.00 | 0.26 | 0.07 | -0.23 | 0.36
86 | 0.19 | 0.04 | -0.02 | 0.10 | 0.15 | 0.04 | 0.26
27 | 0.19 | -0.02 | -0.20 | 0.16 | 0.21 | 0.02 | 0.40
95 | 0.19 | 0.07 | -0.07 | 0.22 | 0.11 | -0.02 | 0.25
16 | 0.18 | -0.01 | -0.12 | 0.10 | 0.19 | 0.09 | 0.29
81 | 0.18 | 0.05 | -0.00 | 0.11 | 0.12 | 0.04 | 0.21
21 | 0.17 | 0.06 | -0.02 | 0.14 | 0.09 | -0.02 | 0.20
74 | 0.16 | 0.09 | -0.03 | 0.20 | 0.07 | -0.16 | 0.29
30 | 0.16 | -0.15 | -0.34 | 0.04 | 0.32 | 0.13 | 0.51
88 | 0.15 | 0.10 | -0.26 | 0.47 | 0.07 | -0.62 | 0.75
15 | 0.14 | 0.03 | -0.09 | 0.16 | 0.11 | -0.07 | 0.29
73 | 0.14 | 0.07 | -0.02 | 0.17 | 0.06 | -0.11 | 0.22
37 | 0.14 | 0.01 | -0.17 | 0.20 | 0.11 | -0.08 | 0.30
58 | 0.14 | 0.08 | -0.05 | 0.21 | 0.06 | -0.10 | 0.22
68 | 0.13 | 0.01 | -0.18 | 0.19 | 0.14 | -0.15 | 0.43
38 | 0.13 | 0.04 | -0.03 | 0.11 | 0.09 | -0.03 | 0.22
75 | 0.12 | 0.04 | -0.19 | 0.27 | 0.07 | -0.30 | 0.44
79 | 0.12 | -0.02 | -0.24 | 0.21 | 0.16 | -0.20 | 0.52
Table 9:  Observed gap, $\widehat{GPG}_{Qd}^{B}$ and $\widehat{GPG}_{Ud}^{B}$
with the confidence interval for the 33 activities with the highest observed
gaps.

##  6 Conclusions

The first contribution of this paper, and the most important from the
theoretical point of view, is the derivation of a novel approach to estimate
the components of the gender pay gap in small areas. This methodology is
consistent with that used by Eurostat (Leythienne and Ronkowski, 2018), it
also generates confidence intervals for the explained and unexplained part of
the GPG, includes a bias correction and considers linear mixed models and a
selection of models option. Moreover, this novel approach is robust against
potential differences in unobserved variables.

The second, and most important contribution from the practical point of view,
is that the application of the methodology to estimate wage discrimination in
economic activities in Galicia reveals important differences in the region and
among activities.

In Galicia, $25\%$ of the pay gap is explained by the differences between
characteristics and $75\%$ is due to discrimination. In other words, if there
were no discrimination and the characteristics of men and women had been paid
at the same prices, the wage differential would be reduced to 2.8%.

Regarding the activities, a high pay gap due to differences in the
characteristics is shown in such activities as 66 and 14, while a high pay gap
due to discrimination is shown in other activities, such as 28, 19, 17, among
others.

The existence of salary differences by sex continues to be a real problem,
according to the results obtained. The $GPG$ indicator, together with the
explained and unexplained gap by economic activities, allow for a better
identification and interpretation of the causes of the gender pay gap. As a
consequence, policy actions can be better targeted.

Finally, to determine whether two workers have the same qualification and
value, useful information, such as the total work experience, knowledge of
languages or the awards or publications obtained, among many others, is still
missing. The inclusion of this information as an explanatory variable would
likely derived in more accurate pay gap estimators.

##  7 Appendix

Tables 10 and 11 show the codes and the description of the Statistical
Classification of Economic Activities in the European Community (NACE09 Rev.2)
at the division level, excluding the primary sector.

NACE09 | Description
---|---
6 | Extraction of crude petroleum and natural gas
7 | Mining of metal ores
8 | Other mining and quarrying
9 | Mining support service activities
10 | Manufacture of food products
11 | Manufacture of beverages
12 | Manufacture of tobacco products
13 | Manufacture of textiles
14 | Manufacture of wearing apparel
15 | Manufacture of leather and related products
16 | Manufacture of wood and of products of wood and cork, except furniture
17 | Manufacture of paper and paper products
18 | Printing and reproduction of recorded media
19 | Manufacture of coke and refined petroleum products
20 | Manufacture of chemicals and chemical products
21 | Manufacture of basic pharmaceutical products and pharmaceutical preparations
22 | Manufacture of rubber and plastic products
23 | Manufacture of other non-metallic mineral products
24 | Manufacture of basic metals
25 | Manufacture of fabricated metal products, except machinery and equipment
26 | Manufacture of computer, electronic and optical products
27 | Manufacture of electrical equipment
28 | Manufacture of machinery and equipment n.e.c.
29 | Manufacture of motor vehicles, trailers and semi-trailers
30 | Manufacture of other transport equipment
31 | Manufacture of furniture
32 | Other manufacturing
33 | Repair and installation of machinery and equipment
35 | Electricity, gas, steam and air conditioning supply
36 | Water collection, treatment and supply
37 | Sewerage
38 | Waste collection, treatment and disposal activities; materials recovery
39 | Remediation activities and other waste management services
41 | Construction of buildings
42 | Civil engineering
43 | Specialised construction activities
45 | Wholesale and retail trade and repair of motor vehicles and motorcycles
46 | Wholesale trade, except of motor vehicles and motorcycles
47 | Retail trade, except of motor vehicles and motorcycles
49 | Land transport and transport via pipelines
50 | Water transport
51 | Air transport
52 | Warehousing and support activities for transportation
53 | Postal and courier activities
55 | Accommodation
56 | Food and beverage service activities
Table 10:  List of codes for the NACE09 at division level, excluding the primary sector. NACE09 | Description
---|---
58 | Publishing activities
59 | Motion picture, video and television programme production, sound recording and music
60 | Programming and broadcasting activities
61 | Telecommunications
62 | Computer programming, consultancy and related activities
63 | Information service activities
64 | Financial service activities, except insurance and pension funding
65 | Insurance, reinsurance and pension funding, except compulsory social security
66 | Activities auxiliary to financial services and insurance activities
68 | Real estate activities
69 | Legal and accounting activities
70 | Activities of head offices; management consultancy activities
71 | Architectural and engineering activities; technical testing and analysis
72 | Scientific research and development
73 | Advertising and market research
74 | Other professional, scientific and technical activities
75 | Veterinary activities
77 | Rental and leasing activities
78 | Employment activities
79 | Travel agency, tour operator and other reservation service and related activities
80 | Security and investigation activities
81 | Services to buildings and landscape activities
82 | Office administrative, office support and other business support activities
84 | Public administration and defence; compulsory social security
85 | Education
86 | Human health activities
87 | Residential care activities
88 | Social work activities without accommodation
90 | Creative, arts and entertainment activities
91 | Libraries, archives, museums and other cultural activities
92 | Gambling and betting activities
93 | Sports activities and amusement and recreation activities
94 | Activities of membership organisations
95 | Repair of computers and personal and household goods
96 | Other personal service activities
97 | Activities of households as employers of domestic personnel
98 | Undifferentiated goods- and services-producing activities of private households for own use
99 | Activities of extraterritorial organisations and bodies
Table 11:  List of codes for the NACE09 at division level, excluding the
primary sector (continuation).

##  8 Bibliography

Akaike, H. (1973). Information theory and the maximum likelihood principle. In
International Symposium on Information Theory, pages 267-281, Budapest.
Akademiai Kiado.

Aláez, R., Ullibarri, M. (2000). Discriminación salarial por sexo: un análisis
del sector privado y sus diferencia regionales en España. ICE, Revista de
Economía, 789, 117-138.

Aláez, R., Longás, J.C., Ullibarri, M. (2001). Visualising gender wage
differences in the European Union. Gender, Work and Organization.

Aláez, R., Longás, J.C., Ullibarri, M. (2003). Diferencias salariales en
España: un análisis sectorial/regional. Investigaciones regionales, 3, 5-24.

Anastasiade, M.C. and Tillé, Y. (2017a). Decomposition of gender wage
inequalities through calibration: Application to the Swiss structure of
earnings survey. Survey Metodology, 43:211-234.

Anastasiade, M.C. and Tillé, Y. (2017b). Gender wage inequalities in
Switzerland: the public versus the private sector. Statistical Methods and
Applications, 26 (2) 293-316.

Battese, G. E., Harter, R. M. and Fuller, W. A. (1988). An error-components
model for prediction of county crop areas using survey and satellite data.
Journal of the American Statistical Association, 83(401),28-36.

Berk, R., Brown, L. Buja, a. Zhang, K. and Zaho L. (2013). Valid post-
selection inference. The Annals of Statistics, 41(2), 802-837.

Blinder, A.S. (1973). Wage Discrimination: Reduced Form and Structural
Estimates. Journal of Human Resources, 8(4), 436-455.

De la Rica, S., Dolado, J.J., Llorens, V. (2008). Ceilings or floors: gender
wage gaps by education in Spain. J. Popul Econ, 21, 751-778

Dempster,A. P., Rubin D. B. and Tsutakawa R.K. (1981). Estimation in
covariance components models. Journal of the American Statistical Association,
76, 341-353.

Fernández, M., Montuenga, V. M., Romeu, A E. (2000). Wage differencials and
non competitive behaviour in the Spanish industry labor market. Documentos de
traballo análise económica, n 9, IDEGA.

Fan, Y., and Li, R. (2012). Variable selection in linear mixed effects models.
Annals of statistics, 40(4), 2043-2068.

Fortin, N., Lemieux, T., and Firpo, S. (2011). Decomposition methods in
economics. In Ashenfelter, O. and Card, D., editors, Handbook of Labor
Economics, volume 4, 1-102. Elsevier.

Graf, M. , Marín, M. and Molina, I. (2019). A generalized mixed model for
skewed distributions applied to small area estimation. TEST, 28, 565-597.

Hlavac, Marek (2018). oaxaca: Blinder-Oaxaca Decomposition in R. R package
version 0.1.4. https://CRAN.R-project.org/package=oaxaca.

Jann B (2008). The Blinder-Oaxaca Decomposition for Linear Regression Models.
Stata Journal, 8(4), 453-479.

Jiang, J. (2007). Linear and Generalized linear Mixed Models and Their
Applications. Springer Series in Statistics.

Juhn, C., Murphy, K. M., and Pierce, B. (1993). Wage inequality and the rise
in returns to skill. Journal of Political Economy, 101(3):410-442.

Lee, J. D., Sun, D. L., Sun, Y. and Taylor J. E. (2016). Exact post-selection
inference with application to the lasso. The Annals of Statistics, 44(3),
907?927

Leythienne, D. abd Ronkowski P. (2018). A decomposition of the unadjusted
gender pay gap using Structure of Earnings Survey data. Population and social
conditions. Statistical Working Papers. Publications Office of the European
Union.

Lombardía, M.J., López-Vizcaíno, E., and Rueda, C. (2017). Mixed generalized
akaike information (xGAIC) for small area models. Journal of Royal Statistical
Society Series A, 180, 1229-1252.

Moral-Arce, I., Sperlich, S., Fernandez-Sáenz, A.I., Roca, M.J. (2011). Trends
in the Gender Pay Gap in Spain: A Semiparametric Analysis. J. Labor Res, 33,
173?195.

Oaxaca R.L. (1973). Male-Female Wage Differentials in Urban Labor
Markets.International Economic Review, 14(3), 693-709.

Oaxaca R.L., Ransom M.R. (1994). On the discrimination and the decomposition
of wage differentials. Journal of Econometrics, 61, 5-21.

Oaxaca R.L., Ransom M.R. (1998). Calculation of approximate variances for wage
decomposition differentials. Journal of Economic and Social Measurement,
81(1), 154-157.

Opsomer, J., Claeskens, G., Ranalli, M., Kauermann, G., and Breidt, F. (2008).
Non-parametric small area estimation using penalized spline regression.
Journal of Royal Statistical Society Series B, 70, 265-286.

Popli, G.K. (2013). Gender wage differentials in Mexico: a distributional
approach. Journal of the Royal Statistical Society Series A, 176(2), 295-319.

Rao, J.N.K. and Molina, I. (2015). Small area estimation, Second Edition. J.
Wiley, Hoboken (NJ).

Rueda, C. and Lombardía, M.J. (2012). Small area semiparametric additive
monotone models. Statistical Modelling, 12, 527-549.

Tibshirani, R.J., Taylor J. Lockhart R. and Tibshirani R. (2016). Exact Post-
Selection Inference for Sequential Regression Procedures. Journal of the
American Statistical Association, 111(514), 600-620.


# Analysis & Computational Complexity Reduction of
Monocular and Stereo Depth Estimation Techniques

Rajeev Patwari
<EMAIL_ADDRESS>The authors contributed equally Varo Ly 11footnotemark:
1
<EMAIL_ADDRESS>
###### Abstract

Accurate depth estimation with lowest compute and energy cost is a crucial
requirement for unmanned and battery operated autonomous systems. Robotic
applications require real time depth estimation for navigation and decision
making under rapidly changing 3D surroundings. A high accuracy algorithm may
provide the best depth estimation but may consume tremendous compute and
energy resources. A general trade-off is to choose less accurate methods for
initial depth estimate and a more accurate yet compute intensive method when
needed. Previous work has shown this trade-off can be improved by developing a
state-of-the-art method [12] to improve stereo depth estimation.

We studied both the monocular and stereo vision depth estimation methods and
investigated methods to reduce computational complexity of these methods. This
was our baseline. Consequently, our experiments show reduction of monocular
depth estimation model size by $\sim$75% reduces accuracy by less than 2%
(SSIM metric). Our experiments with the novel stereo vision method [12] show
that accuracy of depth estimation does not degrade more than 3% (three pixel
error metric) in spite of reduction in model size by 20%. We have shown that
smaller models can indeed perform competitively. The source code is available
at https://github.com/rajeevpatwari/AnyNet and
https://github.com/rajeevpatwari/Monocular-Depth

##  1 Introduction

The current state of the art autonomous driving and real time navigation
decision making tasks require visual scene understanding as an addendum to
radio or other sensing techniques. Depth estimation in computer vision and
robotics is a necessary and crucial step. This is done today using stereo
vision depth techniques [5], monocular depth techniques [14], and hybrid
approaches [2].

##  2 Background and Related Work

In our project, we implemented monocular and stereo depth estimation using
supervised machine learning techniques. Our motivation for the models chosen
were to choose models that could provide adequate accuracy while also allowing
less computation which would enable depth estimation on lower latency and
lower power systems in real world applications.

###  2.1 Monocular Depth

Monocular depth estimation from a single image can be accomplished using
classical supervised learning autoencoder techniques as in [3], [14], and
[10]. The network in [14] was developed to provide both accuracy and low
latency. Improvements to the speed and computational complexity of the these
models were accomplished by utlizing an efficient pretrained encoder such as
MobileNet and pruning the network using NetAdapt. In our studies we utlized a
basic CNN for both the enoder and decoder similar to [3].

###  2.2 Stereo Depth

In stereo depth estimation using supervised machine learning techniques,
images from two cameras and ground truth disparity data are used to train a
neural network to estimate the disparity. In [5] many state of the art methods
were reviewed. The process of stereo depth estimation using these techniques
generally consisted of feature extraction from both images, feature matching
between the images, disparity estimation, and final refinement of the
disparity map.

##  3 Technical Approach

###  3.1 Monocular Depth

####  3.1.1 Monocular Depth Architecture

Figure 1: Monocular Depth Estimation Model

For monocular depth estimation, we reproduced a model based on [3]. This
consisted of creating a neural network as shown in Fig 1 consisting of four
downscaling blocks, a bottleneck block, and 4 upscaling blocks.

####  3.1.2 Monocular Evaluation Metric

We used three loss functions in our monocular depth estimation technique. The
following equation shows L1 loss used for model optimization.
$L1_{loss}=\sum_{i=1}^{n}|d_{i,true}-d_{i,pred}|$ In addition to L1 loss, we
also used depth smoothness and SSIM [13] per Tensorflow. The model training
involved minimizing all three losses.

We evaluated the accuracy of our generated depth maps by comparing them with
the target depth maps using the structural similarity index metric.

###  3.2 Stereo

####  3.2.1 Stereo Model Architecture

Figure 2: AnyNet [12] Model

Based on our review of [5], we chose to reproduce a stereo vision technique
using the AnyNet [12] model. We chose this model as it was fast, accurate, and
one of the computationally less intensive models. Figure 2 depicts the
architecture of the model which consists of the typical components as
previously discussed, except they are combined over 4 different stages.

The U-Net feature extractor is an approximately 20 layer CNN that extracts the
features from both the left and right images at various scales (1/16, 1/8, and
1/4) of the images. The warping block matches the features between the left
and right images. The disparity network CNN for each stage is approximately 7
layers. For the first stage since the image scale is so small, there is no
feature matching before the disparity network CNN. The refinement component is
accomplished by an additional CNN called the spatial propagation network,
SPNet, which is applied to the predicted disparity map following stage 3 to
further improve the results. As part of our studies, we varied the size and
complexity of SPNet and assessed it’s impact on the accuracy of the results.

####  3.2.2 Stereo Evaluation Metric

The smooth L1 loss was during training for stereo depth estimation. The smooth
L1 loss corresponds to the L2 loss when the difference between the outputs is
below a specified parameter and corresponds to the L1 loss when it is above
this chosen value. This can be tuned as a hyper parameter, $\beta$.

$L1_{smooth\\_loss}=$
$\sum_{i=1}^{n}\begin{cases}\frac{1}{2\beta}(d_{i,true}-d_{i,pred})^{2}&\text{if\;}|d_{i,true}-d_{i,pred}|<\beta\\\
|d_{i,true}-d_{i,pred}|-.05\beta&\text{otherwise}\end{cases}$

To assess our stereo model results, we utilized a three pixel error criteria
as described in [12] and supported code
https://gist.github.com/MiaoDX/8d5f49c2ccb39d7f2cb8d4e57c3ab752.

##  4 Experiments and Results

###  4.1 Dataset

####  4.1.1 Monocular Dataset

KITTI [4], NYU-V2 [7], DIODE [11] are three datasets available for monocular
depth estimation, but all three datasets are approximately 100 GB. To reduce
training time, we have used a validation set from the DIODE dataset which
comprises of 325 indoor depth maps and 446 outdoor depth maps. And each data
sample comprises of RGB image, depth map, and depth validity mask. Each image
is 1024x768, however to reduce the training time we reduced the image size to
256x256 and consequently clipped the depth map as well that is provided as a
label to the network. Using this technique each epoch takes approximately 10
minutes to train on the GTX1070 and V100.

In Fig 3, 3 samples are shown from DIODE dataset’s test set. The dataset was
split between training and validation sets. Validation set was used to
evaluate the model’s performance on unseen data.

Figure 3: DIODE dataset. Left:input and Right:Depth map

Fig 4 shows a point cloud representation of an input image using its
corresponding depth map.

Figure 4: 3D point cloud representation of input

####  4.1.2 Stereo Dataset

KITTI [4], NYU-V2 [7], Scene Flow [6] are available datasets for stereo depth
estimation. Each dataset is on the order of 100GB in size. We picked Monkaa, a
subset of the Scene Flow dataset, as the target dataset for our study allowing
us enough time to perform multiple ablation studies and model architecture
evaluation. It contains 9000 synthetic data samples and corresponding
disparity maps. We split the data into train and test sets with a 90-10 split.
The figure 5 shows reference example of a single data sample from the dataset.

Figure 5: Monkaa Stereo dataset. Left, Right input images and Left Disparity
Map

###  4.2 Baseline

Our baseline is to reproduce the existing results as a first step. In the
subsections below, we describe methods used to reproduce both Monocular and
Stereo Vision based depth estimation techniques. We also show some ablation
studies done to improve accuracy of predictions.

####  4.2.1 Monocular Depth

We reproduced the chosen monocular depth estimation model and evaluated
different types of loss functions. Monocular depth estimation was built using
Tensorflow and Keras framework [1]. The reference publication describes using
structural similarity (SSIM), depth smoothness and L1 loss as three metrics.
We performed training and hyper-parameter search for combination of these loss
functions. Fig 6 and Fig 7 below show the convergence over time. We also
varied batch sizes from 16, 32, 64.

We used SSIM as an evaluation metric for measuring similarity between target
and predicted data. Figure 12 shows results.

Figure 6: Monocular Depth - Train Loss Figure 7: Monocular Depth - Validation
Loss Figure 8: Monocular Depth - Input, Target, Predicted

We observed that using SSIM as only loss metric resulted in better model
convergence, as visualized by training and validation loss curves in figures 6
and 7. However, visually, the results were not on part with using all three
losses in optimization. We found this observation to align with our
expectation. SSIM is generally lenient on errors in nearest neighbors,
however, adding depth smoothness and L1 depth loss would provide a tighter
correspondence between targets and predictions.

####  4.2.2 Stereo Depth - AnyNet

We reproduced the AnyNet learning based stereo vision depth estimation model
using Pytorch [8]. Training process is computationally intensive and took 4
hours on GTX1070 and 30min on V100 GPUs. We were also able to use a larger
dataset on this model compared to the Monocular depth model. The dataset was
split 90%-10% for training and test sets randomly. Figure 9 shows training
convergence of 2 variations of the AnyNet model, with and without Spatial
Propagation Network (SPNet). The publication [12] describes that the model
performs better with SPNet. Figures 10 and 9 shows losses at various stages,
described in the model shown in figure 2. The idea is that the user would be
able to truncate the model at stage 1, 2, 3 with a slight loss in accuracy but
a considerably smaller model size.

The publication [12] uses three pixel error for measuring disparity
differences between target and predicted disparities. We have implemented the
same metric to evaluate accuracy of the model.

Fig 11 shows predicted disparity map. Three pixel error accuracy is shown in
the figure 16 and table 3.

Figure 9: AnyNet - Train Loss Figure 10: AnyNet - Test Loss Figure 11: AnyNet
- Left, Right inputs; Target & Predicted disparity

###  4.3 Complexity Reduction

Once we established our baseline with both the models, our subsequent goal was
to evaluate computational complexity of the baseline models and investigate
methods to reduce the model’s computational complexity. Model size and number
of parameters are our metrics to measure reduction in compute usage. On an
embedded device mounted on a battery power autonomous system, reduction in
model size can boost the efficiency of the whole system. Our evaluation metric
remains the same, i.e. lower loss, lower SSIM for Monocular Depth model and
lower three pixel error for Stereo Depth.

The subsections below describe experiments and results on reducing
computational complexity of the models.

####  4.3.1 Monocular Depth

Monocular depth model has an autoencoder structure with 4 Down-scale blocks, 1
bottleneck block and 4 Up-scale blocks. Table 1 describes structure of the
model (Conv2D is a 2-D convolution; BN is batch normalization; Act is
activation function).

Table 1: Monocular Depth - Model structure  Block | Structure
---|---
Down-scale | 2x(Conv2D-BN-Act)+Pool
Bottleneck | 2x(Conv2D-Act)
Up-scale | 2x(Conv2D-BN-Act)

The filter sizes for the monocular model described in table 1 are in increased
order 16-32-64-128-256. The higher order down-scale and upscale blocks are
computationally expensive. Experiments were done in order to reduce the
dependency on higher order filters and evaluate model accuracy, with reference
to the baseline.

The table below 1 describes three model variants, the first row (Model-1) is
the default monocular depth model. The second row (Model-2) is model variant
by removing the higher order down scale and upscale blocks. This variation
reduces the number of model parameters drastically, as shown in the table. In
both Model-1 and Model-2, the activation function used is default
LeakyRelU(0.2). Activation function was changed to Swish based on [9] to
investigate if it improves accuracy. This variant is shown in last row
(Model-2 Swish) of the table 2. The last column of the table shows SSIM
evaluation between target depth map and predicted depth map after training for
20 epochs and model convergence was achieved.

The number of model parameters reduces from 1966467 to 489091, which is a
75.17% reduction in model size!

Table 2: Monocular Depth - Model variants, sizes & accuracy  Model Structure | Activation | Parameters | Test SSIM
---|---|---|---
4-1-4 | LeakyReLU | 1966467 | 0.9895
3-1-3 | LeakyReLU | 489091 | 0.9903
3-1-3 | Swish | 489091 | 0.9871

For the given dataset, reducing the model structure from 4-1-4 to 3-1-3
reduced model parameters by more than two thirds and yet maintained relatively
acceptable performance with respect to baseline.

The figure 12 shows SSIM comparison of 6 random data samples from test set.
The figure 13 shows visual comparison of the depths maps obtained by 3
different models from 2.

Figure 12: Monocular Depth - Models Comparison

The figure 13 shows visual data for the test samples 0, 2, 3, whose
corresponding SSIM is shown in the figure 12.

Figure 13: Monocular Depth - Input, Target, Model-1, Model-2, Model-2 Swish

Based on this, we can conclude that the model complexity can be reduced
without much impact to predicted depth map as shown by Model-1 vs Model-2. An
interesting observation is that replacing LeakyReLU with another advanced
activation function Swish did not improve the accuracy of the model.

Another important observation can be made that although SSIM measured on the 3
models, as shown in 1 and 12 is close to 1 and very close to each other, the
depth maps shown by figure 13 show huge difference in qualitative perception
of depth map. This can explain that SSIM alone cannot be a good loss function
for training the Monocular Depth Estimation model. Thus, we can conclude that
although the models converged acceptably during initial baseline ablation
studies, predicted depth maps did not resemble target.

###  4.4 Stereo Depth - AnyNet

AnyNet model’s structure is mainly comprised of 3 sections as shown in 2.
First section is a UNet for feature detection of left and right images
separately. Second section is a 4D cost volume reduction. Third section is
Spatial Propagation Network (SPN). SPN has larger number of channels in the
filters. As UNet is crucial part of the model to detect disparity features,
this is not primary focus to reduce complexity of the model. Furthermore, [12]
concretely suggests using larger models instead of using UNet if more accuracy
is needed at the expense of larger model. Thus, we started by studying effect
of SPN on the model accuracy.

Similar to exercise done on previous model, our metric of reducing complexity
of AnyNet is to reduce the number of model parameters, which in turn is a
direct metric of reduction in compute and memory resources required.

We compared 5 model variants of AnyNet in Table 3. Each model’s total number
of parameters and average three pixel error across entire test are shown in
the table. The lesser the three pixel error, the higher the accuracy of the
model. Fig 16 shows distribution of three pixel error across all the samples
in test set ( 800 samples). We do observe some outliers but it is promising to
see that all variants of the model have very low error.

Table 3: AnyNet - Model variants, sizes & accuracy Model Structure | Parameters | Three Pixel Error
---|---|---
No SPN | 34629 | 0.2994
SPN 1 channel | 34827 | 0.3048
SPN 2 channels | 35277 | 0.3193
SPN 4 channels | 36933 | 0.3264
SPN 8 channels | 43269 | 0.3178

The figures 14 and 15 show training and test loss convergence for these 5
model variants.

Figure 14: AnyNet - Train Loss for 5 model variants  Figure 15: AnyNet - Test
Loss for 5 model variants  Figure 16: AnyNet - Three Pixel Error on test set
for 5 model variants

Figures 17 and 18 showcase visually the disparity maps obtained by training
the 5 variants of AnyNet models. As observed, the difference is relatively low
on model with 1 channel SPN vs the model with 8 channel SPN, with 20% lesser
model size as shown by the distribution of the 25th and 75th quartiles.

Figure 17: AnyNet - Test Sample 0 - Top Row (Left to Right): Left input, Right
input, Target Disparity, Disparity with No SPN Model; Bottom Row (Left to
Right): Disparity with 1 channel SPN, Disparity with 2 channel SPN, Disparity
with 4 channel SPN, Disparity with 8 channel SPN  Figure 18: AnyNet - Test
Sample 1 - Top Row (Left to Right): Left input, Right input, Target Disparity,
Disparity with No SPN Model; Bottom Row (Left to Right): Disparity with 1
channel SPN, Disparity with 2 channel SPN, Disparity with 4 channel SPN,
Disparity with 8 channel SPN

From the before mentioned experiments we can conclude that addition of SPN
improves accuracy slightly compared to without SPN. As shown in 3 we can
choose 1 channel SPN based AnyNet instead of 8 channel SPN with 20% lesser
model size.

##  5 Conclusion

As an initial step, we reproduced learning based methods to estimate depth
using both, Monocular and Stereo Vision techniques. We performed ablation
studies to tune hyper parameters such as learning rate in order to improve
accuracy. We investigated on using various loss functions on monocular depth
estimation. After establishing baselines in both methods, we experimented with
modifying model architecture by changing the neural network layer sizes. We
also changed activation functions to improve accuracy after reducing network
size.

For stereo depth estimation, we started from the novel state-of-the-art
approach of [12] and studied advantages of reducing the spatial propagation
network in reducing model complexity.

In both the methods, we established that reducing model size is a viable
option to reduce computational complexity, that consequently reduces energy
consumption for deploying these techniques in highly constrained navigation
systems. To further improve quality of experiments, future work can use same
dataset on both monocular and stereo depth techniques for a quantitative and
qualitative comparison of both models. However, stereo technique may always
outperform monocular depth technique in accuracy.

After reducing the model sizes, we also tried using multi-precision casting in
both Pytorch and Tensorflow but were unsuccessful in converting the model from
a 64b/32b to a bfloat16 or float16 model. This could be a reasonable future
work that can be undertaken as reducing the data type from 32b to 16b further
reduces compute and memory needs by 50%. For the most part, this would only
need fine-tuning the model. Similarly, the models can be further reduced
computationally by pruning and quantizing but this would hurt the accuracy
furthermore. Thus, quantizing and pruning may need additional investigative
work (FastDepth [14] research used pruning in their method).

## References

  * [1] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mané, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viégas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, and X. Zheng.  TensorFlow: Large-scale machine learning on heterogeneous systems, 2015\.  Software available from tensorflow.org.
  * [2] J. S. Ashutosh Saxena and A. Y. Ng.  Depth estimation using monocular and stereo cues, 2007.
  * [3] V. Basu.  Keras-monocular depth estimation, 2021.
  * [4] A. Geiger, P. Lenz, C. Stiller, and R. Urtasun.  Vision meets robotics: The kitti dataset.  International Journal of Robotics Research (IJRR), 2013.
  * [5] H. Laga, L. V. Jospin, F. Boussaid, and M. Bennamoun.  A survey on deep learning techniques for stereo-based depth estimation.  IEEE Transactions on Pattern Analysis and Machine Intelligence, page 1–1, 2020.
  * [6] N. Mayer, E. Ilg, P. Häusser, P. Fischer, D. Cremers, A. Dosovitskiy, and T. Brox.  A large dataset to train convolutional networks for disparity, optical flow, and scene flow estimation.  In IEEE International Conference on Computer Vision and Pattern Recognition (CVPR), 2016.  arXiv:1512.02134.
  * [7] P. K. Nathan Silberman, Derek Hoiem and R. Fergus.  Indoor segmentation and support inference from rgbd images.  In ECCV, 2012.
  * [8] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala.  Pytorch: An imperative style, high-performance deep learning library.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems 32, pages 8024–8035. Curran Associates, Inc., 2019.
  * [9] P. Ramachandran, B. Zoph, and Q. V. Le.  Searching for activation functions.  CoRR, abs/1710.05941, 2017.
  * [10] S. T.-V. Thatcher Freeman.  Cs231a - training monocular depth estimation models on a budget, 2021\.
  * [11] I. Vasiljevic, N. Kolkin, S. Zhang, R. Luo, H. Wang, F. Z. Dai, A. F. Daniele, M. Mostajabi, S. Basart, M. R. Walter, and G. Shakhnarovich.  DIODE: A Dense Indoor and Outdoor DEpth Dataset.  CoRR, abs/1908.00463, 2019.
  * [12] Y. Wang, Z. Lai, G. Huang, B. H. Wang, L. van der Maaten, M. Campbell, and K. Q. Weinberger.  Anytime stereo image depth estimation on mobile devices.  CoRR, abs/1810.11408, 2018.
  * [13] Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli.  Image quality assessment: From error visibility to structural similarity.  Trans. Img. Proc., 13(4):600–612, apr 2004.
  * [14] D. Wofk, F. Ma, T.-J. Yang, S. Karaman, and V. Sze.  Fastdepth: Fast monocular depth estimation on embedded systems, 2019.