Global simulation envelopes for diagnostic plots in regression models

David I. Warton

School of Mathematics and Statistics and Evolution & Ecology Research Centre,
UNSW Sydney, NSW 2052, Australia

ORCID: 0000-0001-9441-6645

> Residual plots are often used to interrogate regression model assumptions,
> but interpreting them requires an understanding of how much sampling
> variation to expect when assumptions are satisfied. In this paper, we
> propose constructing global envelopes around data (or around trends fitted
> to data) on residual plots, exploiting recent advances that enable
> construction of global envelopes around functions by simulation. While the
> proposed tools are primarily intended as a graphical aid, they can be
> interpreted as formal tests of model assumptions, which enables the study of
> their properties via simulation experiments. We considered three model
> scenarios – fitting a linear model, generalized linear model or generalized
> linear mixed model – and explored the power of global simulation envelope
> tests constructed around data on quantile-quantile plots, or around trend
> lines on residual vs fits plots or scale-location plots. Global envelope
> tests compared favorably to commonly used tests of assumptions at detecting
> violations of distributional and linearity assumptions. Freely available R
> software (ecostats::plotenvelope) enables application of these tools to any
> fitted model that has methods for the simulate, residuals and predict
> functions.
>
> _Keywords: assumption checking, goodness-of-fit test, Monte Carlo
> simulation, normal quantile plot, parametric bootstrap, residuals vs fits
> plot_

##  1 Introduction

When fitting a regression model, residual plots are typically relied upon as a
diagnostic tool to check assumptions are reasonable. However, a challenge for
interpretation is that residuals will always stray from expected trends to
some extent – we would like to know if deviations from expected patterns are
large relative to sampling variation.

One way to compare observed patterns to expected sampling variation is to
construct a global envelope – a region of the diagnostic plot within which we
expect the observed data to always lie if assumptions are satisfied (at a
chosen level of confidence $1-\alpha$). Rosenkrantz (2000) developed an
approach to construct global confidence bands around a quantile-quantile (or
“QQ”-) plot, under the assumption that observations are independent and
normally distributed, by exploiting Mood joint confidence sets for the mean
and variance. Aldor-Noiman et al. (2013) devised a closely related method, but
used simulation to estimate the desired multiple testing correction needed for
global envelopes.

Many regression models do not assume errors are independent and normally
distributed, as for example when analysing discrete data (Figure 1). In such
cases simulation could be used to construct envelopes, repeatedly simulating
new responses from the fitted model, refitting the model and constructing
residuals. Moral et al. (2017) used this approach to generate pointwise
simulation envelopes around half normal plots. Note however that pointwise
envelopes can be quite misleading, as they do not correct for multiple
testing. The chance that not all observations fall inside a pointwise envelope
can be much higher than the nominal rate, so naïve interpretation of pointwise
envelopes could too often lead to the conclusion that assumptions have been
violated, in situations where they have not been violated. What is needed here
is to construct global simulation envelopes – controlling the chance that
_all_ observations are contained inside the simulation envelope, across
datasets that satisfy model assumptions.

This paper describes a method of constructing global simulation envelopes for
residual plots from any fitted model, using the method of Myllymäki et al.
(2017). Their method was designed as a diagnostic tool for point patterns, but
in principle, it can be used to construct global simulation envelopes around
any function. Their method is used in this paper to construct global
simulation envelopes around three commonly used types of residual plot (Figure
1), as a visual tool to aid interpretation of residual plots. Such envelopes
may also be useful as a teaching tool, to assist students in learning when
departures from expected trends are large relative to sampling variation. A
global envelope could also be used in a formal test of assumptions, and while
this is not the primary intention of the proposed tool, it does enable study
of its properties via power simulation. Simulations presented here suggest
that global envelopes on residual plots have power that is competitive with
formal tests designed to interrogate model assumptions. The proposed global
envelope approach has been implemented on R as the plotenvelope function in
the ecostats package (Warton, 2022), and can in principle be applied to _any_
model object that has methods for the simulate, residuals and predict
functions.

##  2 Global envelopes

Assume we have a diagnostic plot which defines a function $T(r)$ over some set
of values $r\in\mathcal{I}$, which could be an interval, or a discrete set of
values. In constructing a $100(1-\alpha)\%$ level global envelope over
$\mathcal{I}$, we want to find functions $T_{\text{low}}(r)$ and
$T_{\text{upp}}(r)$ such that when the assumed regression model is satisfied

$\mathbb{P}\left[\exists
r\in\mathcal{I}:T(r)\notin(T_{\text{low}}(r),T_{\text{upp}}(r))\right]=\alpha$

###  2.1 Constructing global envelopes via simulation

Myllymäki et al. (2017) proposed a method of constructing a global envelope
for $T(r)$ via simulation – making inferences from a set of simulated
realisations of the function $T_{b}(r)$ for $b=2,\ldots,B$, where we
additionally include $T_{1}(r)=T(r)$, the observed function, in this
functional sample. While Myllymäki et al. (2017) were motivated by the problem
of diagnosing goodness-of-fit in point process models, e.g. using Ripley’s
(1976) $K$-function, the approach could equally well be applied to any type of
function. Here we will define functions using common residual diagnostic
plots, and simulate new residuals by simulating data from the fitted
regression model, refitting the regression model, and recomputing the
residuals (a “parametric bootstrap”, Davison and Hinkley, 1997).

The core idea (Myllymäki et al., 2017) is to construct a global statistic that
measures the maximal distance $u$ of each function $T_{b}(r)$ from the centre
of the $B$ functions $T_{0}(r)$, and construct envelopes such that the
$(1-\alpha)$-quantile of $u$ is only exceeded with probability $\alpha$. For
example, we could compute the maximum absolute difference (MAD) measure, as
compared to the mean function, and construct global envelopes based on
absolute distance from $T_{0}(b)$. Specifically, if we let the mean function
be $T_{0}(r)=\frac{1}{B}\sum_{b=1}^{B}T_{i}(r)$, then for each function
$T_{b}(r)$ we compute
$u_{b}=\max\limits_{r\in\mathcal{I}}|T_{b}(r)-T_{0}(r)|$. We can then
construct a $100(1-\alpha)$% global envelope as

$\left(T_{0}(r)-u_{1-\alpha},T_{0}(r)+u_{1-\alpha}\right)$

where $u_{1-\alpha}$ is the $(1-\alpha)$-quantile of the $u_{b}$, computed as
the $(1-\alpha)B$-smallest value of $u_{b}$. Since $u_{1-\alpha}$ is only
exceeded by a proportion $\alpha$ of the $u_{b}$, only a proportion $\alpha$
of the simulated functions $T_{b}(r)$ will stray outside this envelope.

An issue with using MAD is that $T(r)$ may have non-constant variance – in
fact for most diagnostic plots we consider below, we would expect $T(r)$ to
have a smaller variance for intermediate values of $r$. We will compute the
pointwise sample variance function
$\text{var}(T(r))=\frac{1}{B-1}\sum_{b=1}^{B}\left(T_{b}(r)-T_{0}(r)\right)^{2}$
and measure maximal standardized distance between $T_{b}(r)$ and $T_{0}(r)$,
computed as
$v_{b}=\max\limits_{r\in\mathcal{I}}\frac{|T_{b}(r)-T_{0}(r)|}{\sqrt{\text{var}(T(r))}}$.
Hence we construct a Studentized MAD global envelope (Myllymäki et al., 2017):

$\left(T_{0}(r)-v_{1-\alpha}\sqrt{\text{var}(T(r))},T_{0}(r)+v_{1-\alpha}\sqrt{\text{var}(T(r))}\right)$

where as before $v_{1-\alpha}$ is the $(1-\alpha)$-quantile of the $v_{b}$. Up
to Monte Carlo error, this envelope provides an exact global test of a simple
hypothesis, where the distribution of simulated data has been exactly
specified (Myllymäki et al., 2017). It has been noticed empirically that the
method is often conservative for composite hypotheses, where the distribution
we wish to sample from is a function of unknown parameters (Myllymäki et al.,
2017). When using global simulation envelopes to diagnose regression models,
we find ourselves in this latter situation. Software to compute Studentized
MAD global envelopes is available on R in `GET::global_envelope_test` using
the option type="st".

###  2.2 Diagnostic plots in regression models

Now consider the situation where we have univariate responses stored in a
vector $\mathbf{y}$ of length $n$, with the $i$th observation denoted $y_{i}$
for $i=1,\ldots,n$. We can define vectors of linear predictors
$\hat{\boldsymbol{\eta}}=g(\hat{\boldsymbol{\mu}})$ and residuals
$\boldsymbol{e}$, whose $i$th elements are $\hat{\eta_{i}}$ and $e_{i}$. The
residual $e_{i}$ is some pre-defined function of $y_{i}$ and $\hat{\eta_{i}}$,
an ideal choice being a function such that the distribution of residuals does
not vary as a function of linear predictors when assumptions are satisfied.
There are many different types of diagnostic plot that could be constructed.
Below we describe three commonly used plots, and how we define the interval
$\mathcal{I}$ and the function $T(r)$ for each plot:

QQ-plot:

This is a plot of sorted residuals against theoretical normal quantiles. Let
$e_{(i)}$ be the $i$th residual when placed in rank order. In a normal
quantile or “QQ”-plot we plot $e_{(i)}$ against
$z_{i}=\Phi^{-1}\left(\frac{i-0.5}{n}\right)$ for $i=1,\ldots,n$, where
$\Phi(z)=\mathbb{P}(Z\leq z)$ is the cumulative distribution function of the
standard normal distribution $Z\sim\mathcal{N}(0,1)$. Inverting this formula
for $z_{i}$, we get $i=n\Phi(z_{i})+0.5$.

We define as our function the sorted residual, as a function of its
corresponding normal quantile. Specifically, we let $\mathcal{I}$ be the set
of normal quantiles $\\{z_{1},\ldots,z_{n}\\}$, and $T(r)=e_{(n\Phi(r)+0.5)}$.

Residuals vs fits plot:

We plot residuals $e_{i}$ against linear predictors $\hat{\eta}_{i}$, for
$i=1,\ldots,n$. If there is no violation of the mean model, residuals $e_{i}$
should have no trend, as a function of linear predictors. We can diagnose this
by fitting a smoother to the residuals, as a function of linear predictors:

$e_{i}=f(\hat{\eta}_{i})+\epsilon_{i}$

where $\epsilon_{i}$ are errors centred on zero, and $f(\cdot)$ is estimated
using some non-parametric or semi-parametric regression method. We will
estimate $f(\cdot)$ using thin plate regression splines (Wood, 2003), assuming
the $\epsilon_{i}$ are Gaussian with constant variance, using mgcv::gam on R
with maximum likelihood estimation (Wood, 2011, method="ml").

The function we will use as a diagnostic tool for residual vs fits plots is
the smoother defined above, $T(r)=f(r)$, acting on the interval containing the
observed range of linear predictors,
$\mathcal{I}=\left(\text{min}(\hat{\eta}_{i}),\text{max}(\hat{\eta}_{i})\right)$.

Scale vs location plot:

To look for violation of variance assumptions it is common to plot the
absolute value of residuals $|e_{i}|$ against linear predictors
$\hat{\eta}_{i}$, for $i=1,\ldots,n$. If there is no violation of the mean
model, there should be no trend. We can diagnose this by fitting a smoother to
$|e_{i}|$, as a function of linear predictors:

$|e_{i}|=f_{s}(\hat{\eta}_{i})+\epsilon_{i}$

where as previously $\epsilon_{i}$ are errors centred on zero, and we will
estimate $f_{s}(\cdot)$ using mgcv::gam(…, method="ml").

The function we will use as a diagnostic tool for scale-location plots is the
smoother defined above, $T(r)=f_{s}(r)$, acting on the interval containing the
observed range of linear predictors,
$\mathcal{I}=\left(\text{min}(\hat{\eta}_{i}),\text{max}(\hat{\eta}_{i})\right)$.

###  2.3 Global envelopes for residual plots

We have described how to construct a global envelope via simulation, using the
methods of Myllymäki et al. (2017), given a set of simulated functions
$T_{b}(r)$, $b=1,\ldots B$, and we have defined three different types of
function $T(r)$ based on different residual diagnostic plots that could be
used to query goodness-of-fit of a regression model
$\mathcal{M}(\widehat{\boldsymbol{\theta}})$ that has been fitted to sample
data $\mathbf{y}$. We will construct simulated realisations of these functions
by simulating data $\mathbf{y}_{b}$ from the fitted regression model and
recomputing the relevant residual function $T_{b}(r)$, for each
$b=1,\ldots,B$, as follows:

  * •

Simulate data $\mathbf{y}_{b}$ from fitted regression model,
$\mathcal{M}(\widehat{\boldsymbol{\theta}})$

  * •

Refit the regression model to simulated data, to obtain
$\mathcal{M}(\widehat{\boldsymbol{\theta}}_{b})$

  * •

Compute residuals $\mathbf{e}_{b}$ for the simulated data $\mathbf{y}_{b}$
from the refitted model $\mathcal{M}(\widehat{\boldsymbol{\theta}}_{b})$.

  * •

Hence compute $T_{b}(r)$ from resampled residuals $\mathbf{e}_{b}$ and (if
required) the observed linear predictors $\widehat{\boldsymbol{\eta}}$.

This procedure is sometimes referred to as a parametric bootstrap (Davison and
Hinkley, 1997). Note that this procedure requires us to fit the model
$\mathcal{M}(\boldsymbol{\theta})$ to data $B$ times, which can potentially be
quite computationally intensive, if the model cannot be fitted quickly.

In residual vs fits and scale-location plots, note that linear predictors were
not resampled, but were kept fixed at their observed value
$\widehat{\boldsymbol{\eta}}$ as computed using the regression model fitted to
the observed data $\mathcal{M}(\widehat{\boldsymbol{\theta}})$. Initial
simulations considered resampling the linear predictor – computing
$\widehat{\boldsymbol{\eta}}_{b}$ from
$\mathcal{M}(\widehat{\boldsymbol{\theta}}_{b})$ to use in residual plots –
but in noisy data settings, the variability in
$\widehat{\boldsymbol{\eta}}_{b}$ often suppressed the signal in subsequent
residual plots.

###  2.4 Software implementation in the ecostats package

Software to implement this approach is available from the Comprehensive R
Archive Network (R Core Team, 2022) in the ecostats package, which was written
to accompany an introductory text on regression modelling (Warton, 2022). The
key function is plotenvelope, which by default will produce a residual vs fits
plot and a QQ-plot as discussed in this paper, with accompanying 95% global
envelopes constructed via calls to `GET::global_envelope_test`. The user can
control the method of envelope construction and the number of sets of
simulated residuals used. Default settings construct Studentized MAD global
envelopes (type="st") from $B=199$ sets of residuals. A scale-location plot
can be constructed using the argument which=3 (analogous to behaviour of
plot.lm in base R).

As an example of how to use plotenvelope, Brooks et al. (2017a) used
generalised linear mixed models to study counts of salamander abundance across
23 sites, with four replicates at each site, in order to look for an impact of
mining. Below we fit a Poisson log-linear model to the data with a random
intercept for site (using the glmmTMB package), then check assumptions using
all three of the proposed plots:

      library(glmmTMB)
      data(Salamanders)
      m1 <- glmmTMB(count ~ mined + (1|site),family=poisson, data=Salamanders)
      library(ecostats)
      plotenvelope(m1,which=1:3)

We see for this dataset (Figure 1) signs of overdispersion relative to the
Poisson distribution, with points lying above their simulation envelope in the
quantile plot (Figure 1b), and an increasing trend that exceeds its bounds in
the scale-location plot (Figure 1c), for large fitted values.

The ecostats::plotenvelope function can in principle be applied to _any_ model
object that has methods for the simulate, residuals and predict functions. To
date it has been tested on common tools for fitting linear or generalised
linear models, mixed models (Bates et al., 2015a; Brooks et al., 2017a),
generalised additive models (Wood, 2011), and software for multivariate
responses (Wang et al., 2012).

The method used to construct residuals and linear predictors in plotenvelope
can be controlled by the user via the resFunction and predFunction arguments,
which must be functions that take a fitted model as input and are expected to
return residuals and linear predictors, respectively. These arguments default
to the generic residuals and predict functions, respectively, but if
unavailable for a particular model object then resFunction and predFunction
could be used to define suitable alternatives. If desired, these functions
could also be used to construct other types of plots beyond residual vs fitted
value plots, by setting these functions such that they return values beyond
residuals and/or fitted values. For example, we could use this function to
obtain global envelopes around the expected trend on an added variable plot by
setting predFunction to return the predictor of interest.

Figure 1: Example diagnostic plots with global simulation envelopes,
constructed using ecostats::plotenvelope, for the Salamanders data (Price et
al., 2016) available in the glmmTMB package on R (Brooks et al., 2017a). Three
types of residual plot are presented: (a) a residual vs fitted values plot,
(b) a QQ-plot, (c) a scale-location plot. A Poisson loglinear model with a
random intercept was fitted using glmmTMB, and Pearson residuals computed from
the fit. Global envelopes are the shaded areas, not that in (a) the envelope
is so narrow that it is hardly visible. Note from (b) that many residuals on
the QQ-plot fall outside their simulation envelope (taking unusually large
values), and there is an increasing trend on the scale-location plot, which
strays outside its envelope. Both of these trends are suggestive of
overdispersion.

##  3 Simulations

Although global simulation envelopes are proposed here primarily as a
graphical aid, they can also be used to construct formal tests of assumptions.
Simulations have been undertaken to study the performance of such tests, as
compared to some other commonly available tests of model fit, for situations
where we have $n$ univariate responses $y_{i}$, $i=1,\ldots,n$, assumed to be
related to a single predictor $x_{i}$ via a linear function
$\eta_{i}=\beta_{0}+x_{i}\beta_{1}$, and we fit one of three different types
of model for the response $y_{i}$:

  * (a)

A linear model, $y_{i}\sim\mathcal{N}(\eta_{i},\sigma^{2})$, fitted on R using
the lm function.

  * (b)

A Poisson log-linear model, $y_{i}\sim\mathcal{P}(\exp\\{\eta_{i}\\})$, fitted
on R using the glm function.

  * (c)

A Poisson log-linear model with random intercept term,
$y_{i}\sim\mathcal{P}(\exp\\{\eta_{i}+\epsilon_{i\text{ mod }5}\\})$ where
$\epsilon_{j}\sim\mathcal{N}(0,\omega^{2})$ for $j\in\\{0,1,2,3,4\\}$, fitted
on R using the glmmTMB::glmmTMB function. The normally distributed random
intercept that has been added to the linear predictor can take five different
values, which we cycle through repeatedly as we move through the dataset.

The glmmTMB package (Brooks et al., 2017b), used to fit Model (c), is an
alternative to the widely used lme4 package (Bates et al., 2015b). It is used
here because it tends to provide faster, more stable fits (while also having a
broader range of possible variance structures and family arguments). This
extra functionality is made available via Template Model Builder (Kristensen
et al., 2016), an environment for statistical model-building that makes use of
automatic differentiation and C++ for computational expediency.

Standardized residuals were used to evaluate fits of Model (a) via the
rstandard function, which are normally distributed when assumptions are
satisfied. For other models, we used deviance residuals (García et al., 2012),
the default residual output for glm. In the mixed model setting (c), fitted
values were estimated marginally with respect to the random effect. Note that
in both (b) and (c), residuals are not normally distributed when assumptions
are satisfied, motivating the use of simulation to construct global envelopes
for residual plots.

Models (a-c) were each fitted to data simulated under settings similar to what
was assumed by the corresponding model, but under three different
configurations, corresponding to three different types of assumption
violations. Specifically, in each case the linear predictor used in
simulations had the form

$\eta_{i}=\beta_{0}+4x_{i}+\beta_{2}x_{i}^{2},$

but data were simulated under the following three different configurations:

Assumptions satisfied

Data were simulated from Models (a-c) as described previously, with
$\beta_{0}=-2$, $\beta_{2}=0$ and $\sigma=0.25$ (and for Model (c),
$\omega=1$).

Mixture distribution

A 9:1 mixture of responses was simulated as follows:

  * (a)

$y_{i}\sim
0.9\mathcal{N}(\eta_{i},\sigma^{2})+0.1\mathcal{N}(\eta_{i},(4\sigma)^{2})$

  * (b)

$y_{i}\sim
0.9\mathcal{P}(\exp\\{\eta_{i}\\})+0.1\mathcal{P}(4\exp\\{\eta_{i}\\})$

  * (c)

$y_{i}\sim 0.9\mathcal{P}(\exp\\{\eta_{i}+\epsilon_{i\text{ mod
}5}\\})+0.1\mathcal{P}(4\exp\\{\eta_{i}+\epsilon_{i\text{ mod }5}\\})$

with $\beta_{0}=-2$, $\beta_{2}=0$, $\sigma=0.25$ and $\omega=1$ as before.
Note that each of Models (a-c) correctly specify the functional form of their
respective model for the marginal mean – linearity is preserved in (a), and
log-linearity in (b-c) – but distributional assumptions of each model have
been violated.

Quadratic response

Data were simulated from Models (a-c) with $\beta_{0}=1$ and $\beta_{2}=-4$
(and, as previously, $\sigma=0.25$ and $\omega=1$). This produced (on the
linear predictor scale) a parabolic response whose axis of symmetry was at the
centre of the $x_{i}$, and whose values for the $\eta_{i}$ were centred on
similar values to those produced by other simulation scenarios.

This gives a total of nine simulation scenarios, consisting of three models
(a-c) fitted under each of three assumption violation configurations – when
assumptions are satisfied, when distributional assumptions are violated, and
when the mean model is not correctly specified. Under each scenario, 1000
datasets were simulated at each of four different sample sizes –
$n\in\\{10,20,40,80\\}$. The following tests were compared, in terms of their
rejection rate at the $\alpha=0.05$ level:

  * •

Global simulation envelopes around smoothers on a residual vs fits plot or
scale-location plot as described in the previous section (using Studentized
MAD envelopes). We expected the residual vs fits plot to be useful for
detecting assumption violations in Quadratic response simulations, and the
scale-location plot to be more useful for Mixture distribution simulations.

  * •

Global simulation envelopes around a normal quantile (QQ-) plot or a
probability (PP-) plot (using Studentized MAD envelopes). We expected these to
be most useful in Mixture distribution simulations.

  * •

For Model (a), a Shapiro-Wilks test was computed using the shapiro.test
function on R.

  * •

For Models (b-c), the maximized log-likelihood was used as a goodness-of-fit
test statistic, and compared to its null distribution as estimated from
simulated datasets $\mathbf{y}_{b}$ described previously. Lack-of-fit is
indicated if the observed value
$\log\mathcal{L}(\boldsymbol{\theta},\mathbf{y})$ is small compared to values
computed for data simulated under the assumed model,
$\log\mathcal{L}(\boldsymbol{\theta},\mathbf{y}_{b})$.

In all cases, simulation envelopes were constructed using 999 simulated
datasets (in addition to the observed data), and interpreted as a hypothesis
test by asking whether the observed data (or smoother) fell outside its 95%
global envelope at any stage. Because this whole process was repeated 1000
times, for each sample size and simulation setting, simulations were quite
computationally intensive and were undertaken on a computational cluster (UNSW
Sydney, 2022), requiring several thousand hours in total computation time on
processors with 20-60GB RAM.

Results for these nine simulation scenarios are summarized in Figure 2. Note
firstly that the global envelopes all adequately control Type I error (Figure
2, left), although some methods were conservative. Myllymäki et al. (2017)
also saw conservativeness in their simulations, and noted that it is often
observed empirically in Monte Carlo testing of composite hypotheses.

When looking at power of different methods to detect assumption violations,
the most striking result was that the most effective tools were simulation
envelopes around a QQ-plot or a residual vs fits plot, for mixture
distributions (Figure 2, centre) and quadratic response (Figure 2, right),
respectively. Note that global envelopes on a QQ-plot performed comparably to
a Shapiro-Wilks test at finding violations of normality (Figure 2a, centre)
and better than a deviance test at detecting violations of a Poisson
assumption (Figure 2b-c, centre). When there was a quadratic response to
predictors, a global envelope around a smoother on a residual vs fits plot was
by far the most effective method (Figure 2, right) for all three models.

Note that for non-Gaussian models (Figure 2b-c) there was some interaction
between non-normality and non-linearity, with residual vs fits plots having
some ability to detect mixture distributions and a QQ-plot having some power
under non-linearity. Mean and variance parameters are orthogonal in the linear
model, hence it is not surprising that diagnostics designed to target the mean
trend performed poorly when distributional assumptions were violated, and
diagnostics designed to target distributional assumptions performed poorly
when the mean model was incorrect. In the Poisson model, in contrast, a
missing quadratic term introduces overdispersion, which can be detected on a
QQ-plot (Figure 2b-c, right).

Figure 2: Simulation results when fitting Models (a-c) to data simulated under
the assumed model (left column), from a mixture distribution (centre) or with
a quadratic response to the predictor (right). In (a), the grey area is a 95%
global envelope around the observed Type I error of an exact test, at
significance level of 0.05, in this simulation design. Note in particular that
all global simulation envelopes adequately control Type I error (left), that
envelopes around QQ-plots (green solid line) were most effective at finding
violations of distributional assumptions (center), and envelopes around a
smoother on a residual vs fits plot (green dotted line) was most effective at
finding violations of linearity (right).

As a secondary result, note that the scale on which diagnostic tools are
plotted was important – QQ-plots consistently outperformed PP-plots at
detecting assumption violations. This was expected since assumption violations
occurred in the tails of the distribution, and plotting residuals on the
quantile scale puts the emphasis on the tails, whereas a PP-plot transforms
residuals to fit in the unit interval.

A global envelope around the smoother on a scale-location plot was the least
effective of the proposed methods. As expected, it was useful for identifying
distributional assumption violations (Figure 2, centre), but was noticeably
less effective than a QQ-plot, or a Shapiro-Wilks test of normality (Figure
2a, centre), or a deviance-based test for Poisson data (Figure 2b, centre).

##  4 Discussion

In this paper we showed how recent advances in tools for global simulation
envelopes can be used to aid in the interpretation of common residual plots.
Previously, global simulation envelopes had been proposed for QQ-plots
(Rosenkrantz, 2000; Aldor-Noiman et al., 2013), assuming residuals were
Gaussian and independently and identically distributed (“iid”). However, one
or both of the Gaussian and iid assumptions will not be satisfied for
residuals from many types of regression model, an issue which we resolve here
by using simulation to explore the expected distribution of residuals and
residual functions, by simulating data from the true model and recomputing
residuals. Software is provided (the plotenvelope function in the ecostats
package) to compute global simulation envelopes for various types of residual
plot for a broad class of regression models, and illustrated its effectiveness
by simulation for a few common types of regression model.

The intention is that global envelopes be used as a visual tool to help the
analyst determine whether any departures from expected trends are large
relative to what would be expected if assumptions were true. This could be
especially useful as a teaching tool, for students not experienced with
interpreting residual plots – in fact, the motivation for this work arose for
the author while writing a text for ecologists on data analysis (Warton,
2022). While not the primary intention, the global envelope approach can also
be interpreted as a formal test of assumptions, which was helpful as a way to
study the effectiveness of the tool. It was interesting to see in simulations
that the proposed methods tended to compare favorably to common goodness-of-
fit tests, suggesting that using global envelopes on residual plots is
relatively effective as a way to detect violations of assumptions.

Simulations reinforce the idea that different diagnostic plots detect
different types of assumption violations, by design, and as such it is a good
idea to use multiple diagnostic plots when interrogating multiple assumptions.
Beyond the plots discussed here, one could also look at response to
predictors, as advocated for by Cook and Weisberg (1997). The plotenvelope
function could be used to produce plots against a predictor of particular
interest by setting predFunction to be a function that returned the predictor
in question.

The performance of the global envelope tests depend heavily on the type of
residual used to construct plots, the scale on which residuals are plotted,
and if envelopes are constructed around smoothers, the method of smoother
construction, as discussed in more detail below.

The way residuals are constructed is clearly an important consideration when
constructing any residual plot. In preliminary investigations, we found for
example that when fitting a generalized linear mixed model, residual plots had
difficulty detecting violations of distributional assumptions (Figure 2c,
centre) when using unscaled residuals ($y_{i}-\widehat{\mu}_{i}$) rather than
deviance residuals. A potentially useful alternative way to construct
residuals for non-Gaussian models is randomized quantile residuals (Dunn and
Smyth, 1996), which introduce jittering to discrete data to ensure that
residuals computed from the true model will be marginally standard normal, for
any parametric model. Such residuals can be computed on R for generalized
linear mixed models using the statmod::qres function, and can be computed by
simulation from broader classes of distribution using the DHARMa package
(Hartig, 2022).

The scale on which residuals are plotted can also be important to performance
of global envelopes, as evidenced by the differences in performance between
QQ-plots and PP-plots (Figure 2, centre and right). While PP-plots are often
used as a diagnostic tool (Hartig, 2022, for example), they were consistently
outperformed by QQ-plots. These plots differ only in their scale, with PP-
plots mapping residuals to the unit interval (via $\Phi^{-1}(e_{i})$), then
comparing to standard uniform rather than to standard normal quantiles. But
violations of assumptions typically are expressed in the tails of a
distribution, with the occasional very large (in magnitude) residual. When
converted to the unit interval, such outliers become values that are unusually
close to zero or one, which does not provide a strong visual cue. For example,
in Figure 1, a striking feature is the outlying Pearson residual, with a value
of about 15, which had an expected normal quantile closer to 3. On a PP-plot,
this would be expressed as a probability that was equal to one up to machine
error, as compared to an expected probability of 0.998. Mapping residuals to
the unit interval compresses the regions of primary interest to narrow ranges
of values, making it difficult to diagnose violations of assumptions. For this
reason Dunn and Smyth (1996) advocated mapping quantile residuals to the
normal distribution.

Preliminary investigations also suggested that when constructing a global
envelope around a smoother, performance of envelope tests was sensitive to the
method of smoother construction. When using the default mgcv::gam method,
which minimizes a generalized cross validation (GCV) measure, our smoothers
had quite poor performance at detecting assumption violations. The relatively
good performance seen in Figure 2 was achieved through use of maximum
likelihood estimation to fit smoothers (method="ml"), which is noticeably more
computationally intensive, but known to perform better at estimating smooth
functions, largely eliminating the chance of “severe undersmoothing failures”
seen in GCV techniques (Wood, 2011).

Finally, it is worth reminding the practitioner that significant violations of
assumptions are not necessarily of practical importance. Linear models for
example are known to be quite robust to violations of normality due to the
Central Limit Theorem, especially at large sample sizes (see Lumley et al.,
2002, for example). Large sample sizes however lead to narrower simulation
envelopes and a greater chance of detecting assumption violations, in the very
situation where they are less important. It should always be kept in mind by
the practitioner that global envelopes are graphical tools that show how large
deviations are relative to those that might be expected due to sampling error,
but they do not indicate if deviations are large enough to be of practical
importance. Effect size is best measured by looking at the magnitude of
deviations from expected behaviour (Kim and Cribbie, 2018, for example),
rather than just looking at whether these deviations stray outside a global
envelope.

### Disclosure statement

The author reports there are no competing interests to declare.

### Acknowledgements

The author’s work was supported by an Australian Research Council Discovery
Projects scheme (DP180103543). Thanks to anonymous reviewers for suggestions
that improved this manuscript.

## REFERENCES

  * Aldor-Noiman et al. (2013) Aldor-Noiman, S., Brown, L. D., Buja, A., Rolke, W., and Stine, R. A.  “The Power to See: A New Graphical Test of Normality Statistical Computing and Graphics The Power to See: A New Graphical Test of Normality.” (2013).
  * Bates et al. (2015a) Bates, D., Mächler, M., Bolker, B. M., and Walker, S. C.  “Fitting Linear Mixed-Effects Models Using lme4.”  Journal of Statistical Software, 67:1–48 (2015a).
  * Bates et al. (2015b) —.  “Fitting Linear Mixed-Effects Models Using lme4.”  Journal of Statistical Software, 67:1–48 (2015b).
  * Brooks et al. (2017a) Brooks, M. E., Kristensen, K., van Benthem, K. J., Magnusson, A., Berg, C. W., Nielsen, A., Skaug, H. J., Mächler, M., and Bolker, B. M.  “glmmTMB balances speed and flexibility among packages for zero-inflated generalized linear mixed modeling.”  R Journal, 9:378–400 (2017a).
  * Brooks et al. (2017b) —.  “glmmTMB balances speed and flexibility among packages for zero-inflated generalized linear mixed modeling.”  R Journal, 9:378–400 (2017b).
  * Cook and Weisberg (1997) Cook, R. D. and Weisberg, S.  “Graphics for Assessing the Adequacy of Regression Models.”  Journal of the American Statistical Association, 92:490–499 (1997).
  * Davison and Hinkley (1997) Davison, A. C. and Hinkley, D. V.  Bootstrap methods and their application.  1\. Cambridge university press (1997).
  * Dunn and Smyth (1996) Dunn, P. K. and Smyth, G. K.  “Randomized Quantile Residuals.”  Journal of Computational and Graphical Statistics, 5:236–244 (1996).
  * García et al. (2012) García, M., Víctor, B. ., Yohai, J., Ben, M. G., and Yohai, V. J.  “Quantile–Quantile Plot for Deviance Residuals in the Generalized Linear Model.”  Journal of Computational and Graphical Statistics, 13:36–47 (2012).
  * Hartig (2022) Hartig, F.  DHARMa: Residual Diagnostics for Hierarchical (Multi-Level / Mixed) Regression Models (2022).  R package version 0.4.5.
URL https://CRAN.R-project.org/package=DHARMa

  * Kim and Cribbie (2018) Kim, Y. J. and Cribbie, R. A.  “ANOVA and the variance homogeneity assumption: Exploring a better gatekeeper.”  British Journal of Mathematical and Statistical Psychology, 71:1–12 (2018).
  * Kristensen et al. (2016) Kristensen, K., Nielsen, A., Berg, C. W., Skaug, H., and Bell, B. M.  “TMB: Automatic Differentiation and Laplace Approximation.”  Journal of Statistical Software, 70:1–21 (2016).
  * Lumley et al. (2002) Lumley, T., Diehr, P., Emerson, S., and Chen, L.  “The Importance of the Normality Assumption in Large Public Health Data Sets.”  Annual Review of Public Health, 23:151–169 (2002).
  * Moral et al. (2017) Moral, R. A., Hinde, J., and Demétrio, C. G.  “Half-normal plots and overdispersed models in R: The hnp package.”  Journal of Statistical Software, 81 (2017).
  * Myllymäki et al. (2017) Myllymäki, M., Mrkvička, T., Grabarnik, P., Seijo, H., and Hahn, U.  “Global envelope tests for spatial processes.”  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 79:381–404 (2017).
  * Price et al. (2016) Price, S. J., Muncy, B. L., Bonner, S. J., Drayer, A. N., and Barton, C. D.  “Effects of mountaintop removal mining and valley filling on the occupancy and abundance of stream salamanders.”  Journal of Applied Ecology, 53:459–468 (2016).
  * R Core Team (2022) R Core Team.  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria (2022).
URL https://www.R-project.org/

  * Ripley (1976) Ripley, B. D.  “The second-order analysis of stationary point processes.”  Journal of Applied Probability, 13:255–266 (1976).
  * Rosenkrantz (2000) Rosenkrantz, W. A.  “Confidence Bands for Quantile Functions: A Parametric and Graphic Alternative for Testing Goodness of Fit.”  American Statistician, 54:185–190 (2000).
  * UNSW Sydney (2022) UNSW Sydney.  Katana.  PVC (Research Infrastructure), UNSW Sydney (2022).
URL https://doi.org/10.26190/669x-a286

  * Wang et al. (2012) Wang, Y., Naumann, U., Wright, S. T., and Warton, D. I.  “mvabund – an R package for model-based analysis of multivariate abundance data.”  Methods in Ecology and Evolution, 3:471–474 (2012).
  * Warton (2022) Warton, D. I.  Eco-Stats: Data Analysis in Ecology, from $t$-tests to multivariate abundances.  Springer Nature (2022).
  * Wood (2003) Wood, S. N.  “Thin plate regression splines.”  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 65:95–114 (2003).
  * Wood (2011) —.  “Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models.”  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 73:3–36 (2011).


#  FedSS: Federated Learning with Smart Selection of clients

Ammar Tahir  Yongzhou Chen  Prashanti Nilayam

###### Abstract

Federated learning provides the ability to learn over heterogeneous user data
in a distributed manner, while preserving user privacy. However, its current
clients selection technique is a source of bias as it discriminates against
slow clients. For starters, it selects clients that satisfy certain network
and system specific criteria, thus not selecting slow clients. Even when such
clients are included in the training process, they either straggle the
training or are altogether dropped from the round for being too slow. Our
proposed idea looks to find a sweet spot between fast convergence and
heterogeneity by looking at smart clients selection and scheduling techniques.

##  1 Background

Last decade saw an explosion of data driven and machine learning based
applications to solve different problems. This naturally lead to a lot of
fruitful discussion about ownership and control of access to this data. Users
are concerned about privacy of their sensitive data and do not prefer sharing
it with the third party organizations. However, the value of models trained on
this user data is quite high to both the organizations and end users. Thus,
recently we have seen a lot of work on privacy preserving designs of machine
learning models, popularly known as federated learning.

Federated learning is a privacy preserving method of distributed learning over
heterogeneous user data. Federated learning follows the philosophy of
"bringing the code to data, instead of bringing data to code" to address the
above mentioned privacy conerns [7, 19]. The architecture presented by
Bonawitz et. al. [7] consists of a server responsible for selecting client
devices for training from a pool of currently available devices for each
round. Server maintains a copy of global model, which it distributes to the
selected client devices at start of each round. Each client trains the model
on their local data and send the gradients back to the server. The server,
after getting gradients from all the clients calculates the average gradient
using the FedAvg algorithm. This gradient is used to update the global model,
which is then distributed in the next round.

Since gradients can still reveal information about the data, to ensure privacy
further the gradient updates sent by client are preserved with secure
aggregation [8]. Secure aggregation prevents the server from learning
individual clients’ gradients, but the server can learn overall aggregate once
response from all clients is received. This coupled with already synchronous
nature of federated learning, makes the performance of training susceptible to
degradation due to stragglers. Thus, it is desirable from the angle of
performance that the selected devices are almost homogeneous with respect to
the network and compute power. Most of the current designs select clients
randomly from a pool of devices that satisfy certain criteria e.g. minimum 2
GB memory, unmetered network etc [7, 16]. This technique of client selection
introduces explicit bias in the system since factors like device memory and
quality of network are directly linked with socioeconomic status [6, 14, 16].
Thus, it’s critical to improve client selection mechanism to build models that
are void of this explicit bias.

##  2 Related Work

There have been a great emphasis on this problem, and there are solutions
varying from model compression to different strategies of client selection. In
this section, first we discuss different works that highlight and try to
address the problem of bias and then we discuss different current client
selection strategies.

###  2.1 Bias Mitigation Strategies

As we discussed, the existing design of federated learning introduces explicit
bias since parameters based on which devices are selected are linked with
socioeconomic factors[6, 14, 16]. Bias can also come because of non-IID
distribution of data, however this is a deep learning optimization problem and
there are some techniques to circumvent it [28, 17, 10]. To mitigate the
explicit bias discussed above, we need to include slow devices in the training
process. Clearly doing this comes at the cost of slow convergence, owing to
the slow clients straggling the process. However, there have been some work on
reducing compute and network costs for slow client at the expense of accuracy
of the model [26, 17, 15]. For example Li et. al. [17] allows different
devices to perform variable amount of workload depending on their resources,
slow clients can run fewer epochs or use lesser input data. Similarly, Konečnỳ
et. al. [15] compresses model updates in a lossy fashion to reduce
communication costs at the expense of accuracy of the sent parameters. Abay et
al[6] propose a fairness aware regularization in the loss function. However,
we argue that such techniques do not truly address the bias problem because
there is still discrimination in how slow clients are handled. An ideal design
should give equal opportunity to contribute to the global model to all
clients.

###  2.2 Clients Selection Strategies

It’s clear that the heterogeneous communication environments and computation
resources at clients can hamper the overall training speed. To accelerate the
convergence speed under this client heterogeneity, existing work has
investigated to make smarter decisions about clients selection to alleviate
the communication overhead.

In the prevalent random client selection strategy, all clients participate in
the client selection phase in every round. The server uniformly and randomly
selects a subset of clients from the pool[20, 18, 24] for training, and the
chosen clients do multiple iterations of SGD on the local data, given the ML
model and the latest parameters from the server. Finally, the server collects
and aggregates the computed gradients to update the global parameters. Li et.
al. [18] provided a necessary convergence condition for federated learning on
non-iid data with partial clients participation. Ruan et. al. [24] offered a
selection scheme that converges even when devices can flexibly join or leave
the training.

(a) Random Selection (b) FedCS

(c) Ours

Figure 1: (a) The round time in the random selection scheme is determined by
the straggler. (b) While the round time in FedCS is shortest, the client3 and
client4 are always excluded because of their long training time and
transmission delay. (c) Our scheme reaches the best trade off between training
time and data heterogeneity.

Some recent work looks to do client selection based on different criteria like
higher potential for global model convergence or good network conditions. Cho
et. al. [10] revealed that a selection biased towards clients with higher
local loss can increase the speed of convergence. This is because higher loss
indicates higher potential for model improvement. The proposed POWER-OF-CHOICE
algorithm can yield up to 3X faster convergence and 10% higher test accuracy
compared with conventional federated learning with random clients selection.
FedCS[22] requests the resource information, such as wireless network
bandwidth and computational capacities, from selected clients before the
distribution of parameters of the global model. It then only collects the
gradients from clients which can update and upload the parameters within a
deadline. Although these biased clients selection models can facilitate
quicker convergence, it sacrifices the original benefit of federated learning:
the heterogeneity of data.

##  3 Motivation

As evident from the related work, federated learning between clients with
heterogeneous data, device and network can result in prolonged convergence
time. Appropriate clients selection decisions can result in quick convergence.
Convergence time is determined by two factors: i) number of rounds until the
model convergence condition is reached and ii) time duration of each round.
The number of rounds can be reduced by selecting clients that add more value
to the learning e.g. high losses. Whereas duration of each round is determined
by the stragglers, thus selecting devices based on device and network
parameters can reduce convergence time.

To avoid the prolonged round duration time caused by stragglers, we can select
clients with similar training time in each round. Meanwhile, we can randomly
choose clients with different network/computation conditions across rounds to
guarantee data heterogeneity. We understand that slow clients might be in
minority in some cases, where arranging separate round for them can jeopardize
their privacy. To ensure privacy for such rounds we can fill up the round
group with faster clients too. We visualize and compare this strategy with
random selection in Bonawitz et, al. [7] and FedCS [22] in fig1.

##  4 System Design

To reduce bias and also ensure faster model convergence, we propose a smart
approach of clients selection. We first collect the device’s compute captured
by FLOPS (Floating point Operations per Second) along with network conditions
such as uplink and downlink bandwidth [11] of every connected client.

Based on these collected parameters, we categorize clients into different
equal sized clusters, each cluster with clients with similar training time. We
also determine optimal number of clusters, $k$, for a given distribution. $k$
is optimized to minimize training time, while ensuring higher possible degree
of anonymity (high cluster sizes). In every round, the FedSS server chooses a
cluster in a round-robin way and selects clients within it randomly to ensure
equal opportunity of every client.

There are two basic intuitions behind the approach. 1) Reduce bias induced due
to barrier to entry for low bandwidth/ low-end devices, which in turn
minimizes socio-economic bias. By giving fair chance to all the clients, model
gets a better chance to learn from different data distributions. This also
prevents against content farms attacks using uncompromised phones with high
availability and bandwidth. 2) Have more coordinated training rounds with
similar performing clients grouped together. We are less likely to run into a
situation in which overall completion time of a round is longer due to a
fraction of low-performing devices.

###  4.1 Dynamic Client Environment Tracking

From figure 2, the steps in one round of training in FedSS comprise: the smart
clients selection, the distribution of models, the parallel clients training
procedure, the collection of updated models, and the aggregation. Except the
first and the last steps, the time it costs in all other steps is decided by
the clients’ environments. The network conditions, such as bandwidth and
propagation delay, determine the time taken to distribute model and collect
gradients. Meanwhile, the client’s local computing power, reflected by
available CPUs, GPUs and memory, determine the time taken to finish the
training procedure. Therefore, FedSS dynamically tracks and records those
parameters of every connected client, and predicts the approximate time it
takes to accomplish one round.

Figure 2: The workflow to finish one round of synchronous training in FedSS .
The bottleneck includes models distribution, models collection and parallel
clients training, which majorly depend on the client’s computation speed and
network conditions.

Let’s assume the model size is $M$, and total number of floating point
operations in the model is $Flops$. Given the measured network uplink
bandwidth $UL_{i}$, downlink bandwidth $DL_{i}$, client’s FLOPS rate
$FlopsRate_{i}$ and the number of samples at client $S_{i}$ for client
$c_{i}$, the overall round time $T_{i}$ for that client is:

$T_{i}=M/UL_{i}+S_{i}\cdot Flops/FlopsRate_{i}+M/DL_{i}$

There are different faithful methods to measure bandwidth and FLOPS. Measuring
bandwidth has been an active area of research for a long time. There are
active ways to measure bandwidth such as speed tests and probe trains [11].
There are passive mechanisms as well which estimate bandwidth based on network
usage, e.g. client hints [1]. Similarly, there are different benchmarks to
measure FLOPS as a proxy of compute power e.g. Linpack [4] and HPL [3].
Additionally, these estimates of compute power can be improved with real time
observations of time taken to run a training round. The mechanisms to track
mentioned metrics is not the contribution of our work, in fact, we rely on
existing work for this.

###  4.2 Smart Client Selection

After getting the client round time estimate, $T_{i}$ for most of the clients,
we are ready to run our clustering algorithm. The clustering algorithm takes
the number of clusters as an input, and based on it, decide the percentiles of
distribution. For example, if number of clusters is 3, selected percentiles
would be 25, 50 and 75. Clients are sorted into these clusters based on their
$T_{i}$’s mean-squared distance from percentile values. Since the purpose of
this clustering is to group together clients with almost similar training
times, equidistant percentile based centroids serve this purpose reasonably
well.

The clusters constructed this way tend to have uneven distribution of clients.
The clusters with smaller average training time tend to have most number of
clients compared to clusters with higher average training time. This means if
we have a round robin pattern of switching between clusters for each training
round, the clients in larger clusters have smaller probability of selection.
Our earlier solution to address this problem was to construct a weighted round
robin pattern where clusters appeared in the pattern based on their sizes. For
example, if cluster sizes were (50%, 25%, 25%), the weighted pattern would be
(1,1,2,3), where the first cluster is selected for 2 rounds compared to the
other two.

However, based on feedback during midterm report, we realise that having
different sized clusters also means that clients have varying degree of
privacy. We use k-Anonymity metric to argue about degree of privacy.
k-Anonymity is a metric used to measure the degree of anonymity in anonymised
datasets. Given some identifiers (internal and external), one cannot narrow
down some datapoint to fewer than k individuals. In our case, if we have
varrying cluster sizes, clients in different clusters have varying degree of
k-anonymity. Thus, to be fair to all clients, it is desirable to have clusters
of equal size.

We use a simple trick to even out clusters created by our algorithm. As we
observed that the slow client clusters tend to have fewer clients. We pick
slowest client from fast client clusters and add them to the small cluster
until sizes are almost evened out. The inverse of this can also be done if
fast client cluster is smaller, but there we are careful so that shifting
clients from slower cluster does not have a very drastic effect.

This process of evening out has an additional benefit, data heterogeneity
amongst clusters. This reduces the socio-economic data homogeneity that could
arise due to having similar type of clients in same clusters.

###  4.3 Optimal Number of Clusters

As we pointed out in 4.2, standalone client selection algorithm requires the
number of clusters as an argument. However, it can be difficult for operator
to know this. Thus, we solve for optimal value of number of clusters. Having
fewer clusters is desirable from the perspective of privacy, whereas more
clusters results in shorter average round time. However, average round time
has diminishing returns as we keep increasing number of clusters. The optimal
number of clusters would be the point after which we see diminishing returns
on improvement in average round time.

Figure 3: Visualization of the process behind finding the optimal number of
clusters (6 in both cases). Also worth noticing is superior efficiency of
FedSS ’s clustering algorithm compared to KMeans.

Thus, to find optimal number of clusters, we simulate training process given
round time estimates of clients with different value of clusters. This gives
us average round time for the given number of clusters. Then, we use Kneedle
[25] to find the optimal point on the curve of average round time vs.
1/cluster size. The figure 3 shows the curve between average round time and
1/cluster size along with the knee point found by Kneedle algorithm.

Figure 3 also compares the effectiveness of our algorithm as compared to
KMeans. This is simulation of training with 10000 clients and 1000 rounds of
training. We also compared it against DBScan and KDE clustering, which had
almost similar curve as KMeans. The reason our clustering algorithm does
better is because it is optimised to reduce round time. Whereas Kmeans and
other clustering algorithms are specialised to just cluster data, these
clusters may not be optimal with respect to reducing round time. Moreover, it
is difficult for clustering algorithms to construct equal sized clusters.

###  4.4 Overhead

This grouping algorithm has a time complexity of $NlogN$, where $N$ is the
total number of client devices. Finding optimal number of clusters can be
costly if the number of clients is quite high. But this computation can be
bounded by limiting maximum number of clusters to try e.g. for 100 clients, it
is not desirable to have tens of clusters. Moreover, the whole algorithm can
be run in parallel with model training and new schedule can be used from next
round onwards. Network conditions as well as training time can vary over time
depending on factors like competing traffic, memory or compute back-pressure.
Therefore, it is desirable to run grouping algorithm periodically on updated
estimates of $T_{i}$.

In our experimental implementation, we did not implement mechanisms to measure
bandwidth and profile FLOPS because our experiments are simulation based.
However, existing mechanisms to do so are very lightweight and do not have
significant overhead. Similarly, the profiling of FLOPS can be done by
measuring time to train a round.

##  5 Implementation

Figure 4: The RPC-based communication protocol

Based on an open-source federated learning benchmark system leaf[9], we
implemented FedSS by writing a RPC-based communication protocol and scheduling
logic. We emulate the real federated learning by doing distributed
computations in multiple processes on a single machine.

Figure4 illustrates the communication protocol in detail. The coordinator
thread in the server asynchronously sends RPC commands to $K$ selected clients
for training with its local data, while the server thread keeps listening
clients’ requests to upload the updated models. After the client finishes
training, the server thread receives results and increments a shared variable
acked by one. The coordinator thread waits until acked equals $K$ to continue
the aggregation step and finish this round.

To emulate the real variable network conditions in clients, initially we
planned to utilize Mahimahi[21] to emulate a trace-based network. However, the
training time can take 10s of hours to run one baseline using the emulated
network, which is intolerable given the time and experiments resource
limitation. Thus we simulate the network delays based on the internet speeds
data from World Population Review [2].

Since we neither own the required number of mobile phones, nor find an
appropriate mobile-device emulator for experiments, we simulate the training
speed based on the benchmarked FLOPS [5] of top 20 most sold mobile phones in
2020 [27].

Figure 5: The total time taken by FedSS , FedCS and Random Selection to train
2800 rounds. Figure 6: CDF of time taken per round for FedSS , FedCS and
Random Selection.

##  6 Evaluation

We evaluate FedSS against our implementation of FedCS and random client
selection. We train a CNN model on subset of Femnist dataset, distributed in a
Non-IID manner between 20 clients. Initially, we wanted to run our experiments
on a much larger number of clients, but it’s not possible to run so many
instances of concurrent CNN models in our machine (neither given VMs). For
FedSS and random client selection, we set number of clients per round to be 5.
For FedCS, we select 8 clients per round but aggregate updates with first 5
responses. Also note that for system and network heterogeneity, we simulate
delays for different clients. However, to ensure fairness between our system
and benchmarks, same random delays are used for clients across all our
experiments. Similarly, same distribution of dataset between clients is used
across our experiments.

We explored different metrics to evaluate our system. Since our trade-off is
between training time and bias, we primarily focus on metrics to capture these
two.

###  6.1 Training Time

Measuring training time is easy, we can individually measure the time taken to
complete one round of training at the server. This time would include time
taken to send model to all clients concurrently, time taken by the clients to
train and send back the gradients. This process is bottlenecked by the slowest
client accepted by the server for aggregation. Thus, time per each round is
measured by the time taken by the slowest client to receive, train and send
back the model.

Figure 5 shows the total time taken by FedSS and benchmarks to run 2800 rounds
of training. FedCS beats our system by approximately 40%, because in each
round it only waits for the fastest clients to send back the gradients. This
can be validated by looking at the CDF of time taken per round by FedCS in
figure 6. On the other hand Random Selection has the longest training time,
since time taken per round in this case would be straggled by the slowest
client. Again this can be validated by looking at CDF for Random Selection in
figure 6, which is shifted to right. In comparison, FedSS which selects
clients with similar training times together is able to run more than 80% of
the rounds much faster saving us about 40 hours in training.

Figure 7: The average accuracy of the global model for FedSS , FedCS and
Random Selection. Figure 8: The average loss of the global model for FedSS ,
FedCS and Random Selection.

###  6.2 Bias

Before measuring bias, let’s take a look at the performance of global model
after 2800 rounds for each of the selection strategies. There are many well
established evaluation measures for classification models including precision,
recall, accuracy, f1-score, etc. For our evaluation, we have selected accuracy
and f1-score as the measures.

####  6.2.1 Accuracy

Accuracy is the measure to determine the correct prediction ratio for the
given dataset. It is calculated as :

$\frac{TP+TN}{TP+TN+FP+FN}$

Figure 7 & 8 respectively show the average accuracy and loss of global model
across all clients. From the first look, it seems if all the strategies have
similar performance with respect to model validity. However, a deeper look at
the model’s metrics show a clearer picture.

Accuracy for 4 Slowest Clients | Accuracy for 4 Fastest Clients
---|---
Client | FedCS |  | Random
---
Selection
FedSS | Client | FedCS |  | Random
---
Selection
FedSS
1 | 68.76 | 80.12 | 79.81 | 1 | 83.28 | 85.92 | 86.80
2 | 72.34 | 76.17 | 76.59 | 2 | 55.97 | 55.97 | 60.37
3 | 56.68 | 67.51 | 71.33 | 3 | 80.07 | 76.56 | 78.90
4 | 67.32 | 71.25 | 70.86 | 4 | 84.88 | 79.42 | 78.13
Average | 66.27 | 73.76 | 74.65 | Average | 76.05 | 74.47 | 76.05
Table 1: A breakdown of accuracy of model across 4 slowest and 4 fastest
clients for each of the selection strategies.

Table 1 shows the breakdown of accuracy of the model across different kinds of
clients. We notice that while FedCS has an average accuracy as good as Random
Selection and FedSS , it has lower accuracy for slower clients. On the other
hand, on faster clients, its accuracy is higher than both Random Selection and
FedSS . This difference of performance on slow clients captures the bias in
FedCS’s client selection strategy. This difference arises from the fact that
faster clients end up getting a lot more say in the training than the slow
clients.

The 1.5% difference of accuracy for FedSS , although negligible, can be
explained from data heterogeneity point of view. Recall that Random Selection
by design does not have any explicit bias due to device heterogeneity. This is
because Random Selection selects clients randomly and then waits for all of
them to respond before aggregating. Thus, the 1% difference in its accuracy
when compared on faster vs slower clients, can solely be attributed to data
heterogeneity. The fact that FedSS does as good as Random Selection while
reducing the training time by a huge margin shows promise that FedSS can
reduce bias and amortize the cost of slow nodes.

####  6.2.2 F1 Score

F score is a measure to determine model performance by using both, precision
and recall values. Precision is the measure to determine the percentage of
correct results out of all the results for a given class. It is given as :

$TP/(TP+FP)$

Recall is the measure to determine the percentage of correct results out of
all the true results for a given class. It is given as :

$TP/(TP+FN)$

F1 score, also known as balanced F score, is the harmonic mean of precision
and recall.

$2*\frac{precision*recall}{precision+recall}$

For a multi-class classification, F1 score is given as the average of F1 score
of each class based on averaging scheme. For our experiments, we have selected
weighted averaging to deal with class imbalance issue due to non-iid
distribution of data across multiple clients.

F1 score for 4 Slowest Clients | F1 score for 4 Fastest Clients
---|---
Client | FedCS |  | Random
---
Selection
FedSS | Client | FedCS |  | Random
---
Selection
FedSS
1 | 67.10 | 78.08 | 77.22 | 1 | 81.48 | 83.72 | 84.45
2 | 69.00 | 72.58 | 73.51 | 2 | 56.34 | 56.86 | 60.85
3 | 54.04 | 66.32 | 69.97 | 3 | 77.50 | 73.60 | 75.64
4 | 62.28 | 68.02 | 67.69 | 4 | 85.13 | 79.39 | 78.14
Average | 62.60 | 71.25 | 72.09 | Average | 75.11 | 73.39 | 74.77
Table 2: A breakdown of F1 score of model across 4 slowest and 4 fastest
clients for each of the selection strategies. Figure 9: The average F1 score
of the global model for FedSS , FedCS and Random Selection.

Table 2 shows the breakdown of f1-score for the slowest and the fastest
clients trained during our evaluation. The results are consistent with the
accuracy scores. FedSS achieves results similar to FedCS for the fastest
clients and FedSS improves the results by 10% for the slowest clients as
compared to FedCS. Similar to accuracy, the results from FedSS matches to
random selection while taking considerably lesser amount of time to train.

##  7 Future work

We discuss the limitations of FedSS and the future work in this section. While
FedSS achieves the best trade-off between training time and data heterogeneity
in our testbed, we have some restrictions in experiments which need more work
and real world experimentation.

Experimental Setup:  Currently, we are bound to work on our local setup, which
restricts us to be able to run experiments with only 20 clients. Running with
more clients results in system crashing. We believe better performance on
larger number of clients as this improves clustering. Also, it would be really
valuable to be able to run our setup with some real devices, along with
implementations to track client environment.

Handling Churn:  Currently FedSS assumes that all clients join in the
beginning and maintain available during the whole training procedure, which
may not be true in real-world federated learning. Recalling that it takes a
few rounds for FedSS to measure and profile the new client’s network condition
and computation capability, it can be hard for FedSS to track this information
accurately in the scenario where clients join and leave frequently.

To mitigate the possible performance degradation due to the measurement and
profiling delay of new connected clients, we propose to collect the hardware
specifications of clients, such as processor types and memory size, to
initialize the grouping based on some regression models. The rationale behind
it is that the clients with similar hardware should have similar computation
capability under most circumstances (unless many other background tasks
compete for resources). Whenever a new client joins, FedSS takes one RTT to
collect this information and determine its initial group based on the model,
which is still better than random initialization.

Data Heterogeneity:  FedSS treats all clients equally in every selection,
which may not be the optimal solution. Although the selection is unbiased, the
data heterogeneity itself among clients is another important source of
bias[14, 28, 13]. If a subset of clients with similar kind of data are picked
most of the times, model will converge faster but can have severe bias to the
data distribution of that subset and may not be able to capture the true
global data distribution.

To overcome this issue, we want to explore the possibility to introduce the
training loss of individual clients into client selection process.
Specifically in the every round, if the training loss for a particular client
is non-significant, we exclude it from selection for next few rounds. This
will improve the probability of other clients contributing towards the global
model and also prevent model askew.

Scalability and Fault-tolerance FedSS applies only one server in the system,
which induces communication bottleneck and a single point of failure. We plan
to leverage locality-aware multi-server architecture to enhance the
scalability, fault-tolerance and even training performance in the future.

Firstly, a locality-aware server tends to schedule and disseminate its model
to local clients with better network service, which alleviates the
communication bottleneck at the server side and reduces the ratio of
stragglers resulting from network. Secondly, multiple servers can apply a
consensus algorithm to replicate the training state of each other to tolerate
failure. Lastly, for geographically distributed applications such as next work
prediction[12] and geo-local language translation[23], a locality-aware server
can captures and preserves those geo-local characteristics better.

##  8 Conclusion

Federated learning due to its distributed nature of training machine learning
model suffers from the issues of data and device heterogeneity. When we talk
about device heterogeneity, there exists a trade-off between small training
time and bias, as existing schemes end up dropping slower clients. We show
this trade-off by comparing existing mechanisms of client selection. We then
argue that to eliminate bias, it is necessary to make slow clients part of
training. We present FedSS which finds a sweet spot between small training
time and handling device heterogeneity by performing smart selection of
clients.

## References

  * [1] https://developers.google.com/web/funda- mentals/performance/opt-imizing-content-effici- ency/client-hints.
  * [2] https://worldpopulationreview.com/country-rankings/internet-speeds-by-country.
  * [3] http://www.netlib.org/benchmark/hpl/.
  * [4] http://www.netlib.org/linpack/.
  * [5] Passmark software - cpu benchmarks https://www.cpubenchmark.net/.
  * [6] Annie Abay, Yi Zhou, Nathalie Baracaldo, Shashank Rajamoni, Ebube Chuba, and Heiko Ludwig.  Mitigating bias in federated learning.  arXiv preprint arXiv:2012.02447, 2020.
  * [7] K. A. Bonawitz, Hubert Eichner, Wolfgang Grieskamp, Dzmitry Huba, Alex Ingerman, Vladimir Ivanov, Chloé M Kiddon, Jakub Konečný, Stefano Mazzocchi, Brendan McMahan, Timon Van Overveldt, David Petrou, Daniel Ramage, and Jason Roselander.  Towards federated learning at scale: System design.  In SysML 2019, 2019.  To appear.
  * [8] Keith Bonawitz, Vladimir Ivanov, Ben Kreuter, Antonio Marcedone, H Brendan McMahan, Sarvar Patel, Daniel Ramage, Aaron Segal, and Karn Seth.  Practical secure aggregation for federated learning on user-held data.  arXiv preprint arXiv:1611.04482, 2016.
  * [9] S. Caldas, Peter Wu, Tian Li, Jakub Konecný, H. McMahan, Virginia Smith, and Ameet Talwalkar.  Leaf: A benchmark for federated settings.  ArXiv, abs/1812.01097, 2018.
  * [10] Yae Jee Cho, Jianyu Wang, and Gauri Joshi.  Client Selection in Federated Learning: Convergence Analysis and Power-of-Choice Selection Strategies.  arXiv e-prints, October 2020.
  * [11] Bong Dae Choi, Doo Il Choi, Yoonju Lee, and Dan Keun Sung.  Priority queueing system with fixed-length packet-train arrivals.  IEE Proceedings-Communications, 145(5):331–336, 1998.
  * [12] Andrew Hard, Kanishka Rao, Rajiv Mathews, Swaroop Ramaswamy, Françoise Beaufays, Sean Augenstein, Hubert Eichner, Chloé Kiddon, and Daniel Ramage.  Federated learning for mobile keyboard prediction.  arXiv preprint arXiv:1811.03604, 2018.
  * [13] Eunjeong Jeong, Seungeun Oh, Hyesung Kim, Jihong Park, Mehdi Bennis, and Seong-Lyun Kim.  Communication-efficient on-device machine learning: Federated distillation and augmentation under non-iid private data.  arXiv preprint arXiv:1811.11479, 2018.
  * [14] Peter Kairouz, H Brendan McMahan, Brendan Avent, Aurélien Bellet, Mehdi Bennis, Arjun Nitin Bhagoji, Keith Bonawitz, Zachary Charles, Graham Cormode, Rachel Cummings, et al.  Advances and open problems in federated learning.  arXiv preprint arXiv:1912.04977, 2019.
  * [15] Jakub Konečnỳ, H Brendan McMahan, Felix X Yu, Peter Richtárik, Ananda Theertha Suresh, and Dave Bacon.  Federated learning: Strategies for improving communication efficiency.  arXiv preprint arXiv:1610.05492, 2016.
  * [16] T. Li, A. K. Sahu, A. Talwalkar, and V. Smith.  Federated learning: Challenges, methods, and future directions.  IEEE Signal Processing Magazine, 37(3):50–60, 2020.
  * [17] Tian Li, Anit Kumar Sahu, Manzil Zaheer, Maziar Sanjabi, Ameet Talwalkar, and Virginia Smith.  Federated optimization in heterogeneous networks.  arXiv preprint arXiv:1812.06127, 2018.
  * [18] Xiang Li, Kaixuan Huang, Wenhao Yang, Shusen Wang, and Zhihua Zhang.  On the Convergence of FedAvg on Non-IID Data.  arXiv e-prints, July 2019.
  * [19] Brendan McMahan and Daniel Ramage.  Federated learning: Collaborative machine learning without centralized training data, Apr 2017.
  * [20] H. Brendan McMahan, Eider Moore, Daniel Ramage, Seth Hampson, and Blaise Agüera y Arcas.  Communication-efficient learning of deep networks from decentralized data, 2017.
  * [21] Ravi Netravali, Anirudh Sivaraman, Somak Das, Ameesh Goyal, Keith Winstein, James Mickens, and Hari Balakrishnan.  Mahimahi: Accurate record-and-replay for HTTP.  In 2015 USENIX Annual Technical Conference (USENIX ATC 15), pages 417–429, Santa Clara, CA, 2015.
  * [22] Takayuki Nishio and Ryo Yonetani.  Client selection for federated learning with heterogeneous resources in mobile edge.  ICC 2019 - 2019 IEEE International Conference on Communications (ICC), May 2019.
  * [23] Christiane Nord.  Text analysis in translation: Theory, methodology, and didactic application of a model for translation-oriented text analysis.  Number 94. Rodopi, 2005.
  * [24] Yichen Ruan, Xiaoxi Zhang, Shu-Che Liang, and Carlee Joe-Wong.  Towards Flexible Device Participation in Federated Learning for Non-IID Data.  arXiv e-prints, June 2020.
  * [25] Ville Satopaa, Jeannie Albrecht, David Irwin, and Barath Raghavan.  Finding a" kneedle" in a haystack: Detecting knee points in system behavior.  In 2011 31st international conference on distributed computing systems workshops, pages 166–171. IEEE, 2011.
  * [26] Zirui Xu, Zhao Yang, Jinjun Xiong, Jianlei Yang, and Xiang Chen.  Elfish: Resource-aware federated learning on heterogeneous edge devices.  arXiv preprint arXiv:1912.01684, 2019.
  * [27] Yordan.  Top 20 most popular phones in 2020, Dec 2020.
  * [28] Yue Zhao, Meng Li, Liangzhen Lai, Naveen Suda, Damon Civin, and Vikas Chandra.  Federated learning with non-iid data.  arXiv preprint arXiv:1806.00582, 2018.