A model robust sub-sampling approach for
Generalised Linear Models in Big data settings
Amalan Mahendran School of Mathematical Sciences and Centre for Data Science,
Queensland University of Technology

Helen Thompson

School of Mathematical Sciences

and Centre for Data Science,

Queensland University of Technology

James M. McGree

School of Mathematical Sciences

and Centre for Data Science,

Queensland University of Technology

In today’s modern era of Big data, computationally efficient and scalable
methods are needed to support timely insights and informed decision making.
One such method is sub-sampling, where a subset of the Big data is analysed
and used as the basis for inference rather than considering the whole data
set. A key question when applying sub-sampling approaches is how to select an
informative subset based on the questions being asked of the data. A recent
approach for this has been proposed based on determining sub-sampling
probabilities for each data point, but a limitation of this approach is that
appropriate sub-sampling probabilities rely on an assumed model for the Big
data. In this article, to overcome this limitation, we propose a model robust
approach where a set of models is considered, and the sub-sampling
probabilities are evaluated based on the weighted average of probabilities
that would be obtained if each model was considered singularly. Theoretical
support for such an approach is provided. Our model robust sub-sampling
approach is applied in a simulation study and in two real world applications
where performance is compared to current sub-sampling practices. The results
show that our model robust approach outperforms alternative approaches.

Keywords: Large data set, model averaging, optimal sub-sampling

<EMAIL_ADDRESS>
##  1 Introduction

Big data analysis has gained interest in recent years through providing new
insights and unlocking hidden knowledge in different fields of study [18]
including medicine [31], fraud detection [36], the oil and gas industry [29]
and astronomy [48]. However, the analysis of Big data can be challenging for
traditional statistical methods and standard computing environments [37]. [26]
discuss storage and modelling solutions when handling such a large amount of
data. In general, modelling solutions can be grouped into three broad
categories: 1) _sub-sampling methods_ , where the analysis is performed on an
informative sub-sample obtained from the Big data [19, 24, 25, 14, 40, 39, 44,
2, 9, 21, 3, 45]; 2) _divide and recombine methods_ , where the Big data set
is divided into smaller blocks, and then the intended statistical analysis is
performed on each block and subsequently recombined for inference [23, 16, 11,
8, 22]; 3) _online updating of streamed data_ , where statistical inference is
updated as new data arrive sequentially [32, 41]. In recent years, compared to
divide and recombine methods, sub-sampling has been applied to a variety of
regression problems, while online updating is typically only used for
streaming data. In addition, in cases where a large data set is not needed to
answer a specific question with sufficient confidence, sub-sampling seems
preferable as the analysis of the data can often be undertaken with standard
methods. Moreover, the computational efficiency of sub-sampling over analysing
large data sets has been observed for parameter estimation in linear [39] and
logistic [40] regression models. For these reasons, we focus on sub-sampling
methods in this article.

The key challenge for sub-sampling methods is how to obtain an informative
sub-sample that can be used to efficiently answer specific analysis questions
and provide results that align with the analysis of the whole Big data set.
Two approaches for this exist in the literature: 1) randomly sample from the
Big data with sub-sampling probabilities that are found based on a specific
statistical model and objective (e.g., prediction, parameter estimation) [40,
44, 39, 2, 9, 3, 21, 45]; 2) select sub-samples based on an experimental
design [14, 12]. Randomly sampling with certain probabilities (based upon the
definitions of $A$\- or $L$\- optimality criteria, see [5]) is the focus of
this article, and has been applied for parameter estimation in a wide range of
regression problems including softmax [44] and quantile regression [2], and
generalised linear models [40, 39, 9, 3, 45]. In contrast, the approach based
on an experimental design has only been applied for: 1) parameter estimation
in logistic fixed and mixed effects regression models [14]; 2) parameter
estimation and prediction accuracy in linear and logistic regression models
[12].

A key feature of both of the current sub-sampling approaches is that they rely
on a statistical model that is assumed to appropriately describe the Big data.
Given this is a potentially limiting assumption, [47] proposed to select the
best candidate model from a pool of models based on the Bayesian Information
Criterion (BIC) [33]. This was applied to linear models, and resulted in sub-
sampling probabilities that were more appropriate than those based on
considering a single model. Similarly, for linear regression, [34, 27]
explored using space filling and orthogonal Latin hypercube experimental
design techniques to allow for potential model misspecification. In this
paper, we propose that, instead of selecting a single best candidate model for
the Big data, a set of models is considered and a model averaging approach is
used for determining the sub-sampling probabilities. Through adopting such an
approach, it is thought that the analysis goal (e.g., efficient parameter
estimation) should be achieved regardless of the preferred model for the data.
To implement this model robust approach, we consider sub-sampling based on
$A$\- and $L$\- optimality criteria within the Generalised Linear Modelling
framework, and provide theoretical support for using a model averaged approach
based on each of these criteria. Given we consider Generalised Linear Models
(GLMs), our approach should be generally applicable across many areas of
science and technology where a variety of data types are observed. This is
demonstrated through applying our proposed methods within a simulation study,
and for the analysis of two real-world Big data problems.

The remainder of the article is structured as follows. Section 2 introduces
GLMs and the existing probability-based sub-sampling approach of [3] and [45].
Our proposed model robust sub-sampling approach is introduced in Section 3,
which is embedded with a GLM framework. A simulation study is then used to
assess the performance of our model robust approach in Section 4, and two real
world applications are presented. Section 5 concludes the article with a
discussion of the results and some suggestions for future research.

##  2 Background

There are variety of ways Big data can be sub-sampled. In this section, we
focus on the approach where sub-sampling probabilities are determined for each
data point, and the Big data is sub-sampled (at random) based on these
probabilities. Such an approach was first proposed by [40] for logistic
regression problems in Big data settings, and has been extended to a wide
range of regression problems (e.g., [44, 39, 2, 9, 3, 45]). In this section,
we describe such a sub-sampling approach as applied to GLMs based on the work
of [3].

###  2.1 Generalised Linear Models

Let a Big data set be denoted as $F_{N}=(\bm{X}_{0},\bm{y})$, where
$\bm{X}_{0}=({\bm{x}_{0}}_{1},\ldots,{\bm{x}_{0}}_{N})^{T}\in R^{N\times p}$
represents a data matrix based on the Big data set with $p$ covariates,
$\bm{y}=(y_{1},\ldots,y_{N})^{T}$ represents the response vector and $N$ is
the total number of data points. To fit a GLM, consider the model matrix
$\bm{X}=h(\bm{X}_{0})\in R^{N\times(p+q)}$ where $h(.)$ is some function of
$\bm{X}_{0}$ which creates an additional $q$ columns representing, for
example, an intercept and/or higher-order terms. A GLM can then be defined via
three components: 1) distribution of response $\bm{y}$, which is from the
exponential family (e.g., Normal, Binomial or Poisson); 2) linear predictor
$\bm{\eta}=\bm{X}\bm{\theta}$, where
$\bm{\theta}=(\theta_{1},\ldots,\theta_{p+q})^{T}$ is the parameter vector;
and 3) link function $g(.)$, which links the mean of the response to the
linear predictor [28]. Throughout this article, the inverse link function
$g^{-1}(.)$ is denoted by $u(.)$.

A common exponential form for the probability density or mass function of $y$
can be written as:

$f(y;\omega,\gamma)=\exp{\Big{(}\frac{y\omega-\psi(\omega)}{a(\gamma)}+b(y,\gamma)\Big{)}},$
(1)

where $\psi(.),a(.)$ and $b(.)$ are some functions, $\omega$ is known as the
natural parameter and $\gamma$ the dispersion parameter. Based on Equation (1)
the link function $g(.)$ can then be defined as $g(\bm{\mu})=\bm{\eta}$, where
$\mu=E[y|\bm{X},\bm{\theta}]=\mbox{d}\psi(\omega)/\mbox{d}\omega$. A general
linear model, or linear regression model, is a special case of a GLM where
$\bm{y}\sim N(\bm{\mu},\bm{\Sigma})$ and $g(.)$ is the identity link function
$g(\bm{\mu})=\bm{\mu}$, such that $\bm{\mu}=\bm{X}\bm{\theta}$. For logistic
regression, $\bm{y}\sim\mbox{Bin}(n,\bm{\pi})$ and $g(.)$ is the logit link
function $g(\bm{\pi})=\log(\bm{\pi}/(1-\bm{\pi}))$ such that
$g(\bm{\pi})=\bm{X}\bm{\theta}$. Similarly, for Poisson regression,
$\bm{y}\sim\mbox{Poisson}(\bm{\lambda})$ and $g(.)$ is the log link function
$g(\bm{\lambda})=\log(\bm{\lambda})$ such that
$g(\bm{\lambda})=\bm{X}\bm{\theta}$. Note that the dispersion parameter
$\gamma=1$ for the logistic and Poisson regression models.

###  2.2 A general sub-sampling algorithm for GLMs

As described by [3], consider a general sub-sampling approach to estimate
parameters $\bm{\theta}$ through a weighted log-likelihood function for GLMs
(weights are the inverse of the sub-sampling probabilities). A weighted
likelihood function is considered, since an unweighted likelihood leads to
biased estimates of model parameters for logistic regression, see [38]. Define
$\phi_{i}$ as the probability that row $i$ of $F_{N}$ is randomly selected,
for $i=1,\ldots,N$, where $\sum_{i=1}^{N}\phi_{i}=1$ and $\phi_{i}\in(0,1)$. A
sub-sample $S$ of size $r$ is then drawn with replacement from $F_{N}$ based
on $\bm{\phi}=(\phi_{1},\ldots,\phi_{N})$. The selected responses, covariates
and sub-sampling probabilities are then used to estimate the model parameters.
Pseudo-code for this general sub-sampling approach is provided in Algorithm 1.

1 Sampling: Assign $\phi_{i},i=1,...,N,$ for $F_{N}$.

Based on $\bm{\phi}$, draw an $r$ size sub-sample with replacement from
$F_{N}$ to yield
$S=\\{\bm{x}_{l}^{*},y_{l}^{*},\phi_{l}^{*}\\}_{l=1}^{r}=(\bm{X}^{*},\bm{y}^{*},\bm{\phi}^{*})$.

2 Estimation: Based on $S$, find:
$\displaystyle\tilde{\bm{\theta}}=\operatorname*{arg\,max}_{\bm{\theta}}~{}\log{L(\bm{\theta}|\bm{X}^{*},\bm{y}^{*},\bm{\phi}^{*})}\equiv\operatorname*{arg\,max}_{\bm{\theta}}~{}\frac{1}{r}\sum_{l=1}^{r}\frac{y^{*}_{l}u(\bm{\theta}^{T}\bm{x}^{*}_{l})-\psi(u(\bm{\theta}^{T}\bm{x}^{*}_{l}))}{\phi^{*}_{l}}.$

Output: $\tilde{\bm{\theta}}$ and $S$.

Algorithm 1 General sub-sampling algorithm [3]

From Algorithm 1, the first step is to assign sub-sampling probabilities
$\phi_{i}$ to the rows of $F_{N}$. The simplest approach is to assign each
data point an equal probability of being selected. These probabilities could
also depend on the composition of $\bm{y}$, e.g., for binary data, one could
sample proportional to the inverse of the number of successes and failures.
Based on these probabilities, a sub-sample $S$ of size $r$ is then drawn
completely at random (with replacement) from $F_{N}$ to yield
$\\{\bm{x}^{*}_{l},y^{*}_{l},\phi^{*}_{l}\\}_{l=1}^{r}$. Based on this sub-
sample, model parameters $\tilde{\bm{\theta}}$ are estimated via maximising a
weighted log-likelihood function. The estimates $\tilde{\bm{\theta}}$ can then
be considered as estimates of what would be obtained if the whole Big data set
were analysed.

The asymptotic properties of $\tilde{\bm{\theta}}$ based on the general sub-
sampling approach given in Algorithm 1 were derived by [3]. These properties
are outlined below as they form the basis for our extensions to model robust
sub-sampling detailed in Section 3.

###  2.3 Asymptotic properties of $\tilde{\bm{\theta}}$ from the general sub-
sampling algorithm

To explore the asymptotic properties of the estimates obtained from the
general sub-sampling algorithm, [3] made a number of assumptions which are
outlined below. To follow these, note that
$\dot{\psi}(u(\bm{\theta}^{T}\bm{x}_{i}))$ and
$\dot{u}(\bm{\theta}^{T}\bm{x}_{i})$ denote the first-order derivatives of
$\psi(u(\bm{\theta}^{T}\bm{x}_{i}))$ and $u(\bm{\theta}^{T}\bm{x}_{i})$ (with
respect to $\bm{\theta}$), respectively, with two dots similarly denoting the
second-order derivatives. In addition, the Euclidean norm of a vector $\bm{a}$
will be denoted as $||\bm{a}||=(\bm{a}^{T}\bm{a})^{1/2}$.

######  H 1.

Assume that $\bm{X}\bm{\theta}$ lies in the interior of a compact set
$K\in\Omega$ almost surely.

######  H 2.

The regression coefficient $\bm{\theta}$ is an inner point of the compact
domain $\Lambda_{B}=\\{\bm{\theta}\in R^{p}:||\bm{\theta}||\leq B\\}$ for some
constant $B$.

######  H 3.

Central moments condition:
$N^{-1}\sum_{i=1}^{N}|y_{i}-\dot{\psi}(u(\bm{\theta}^{T}\bm{x}_{i}))|^{4}=O_{P}(1)$
for all $\bm{\theta}\in\Lambda_{B}$.

######  H 4.

As $N\rightarrow\infty$, the observed information matrix

$\bm{J_{X}}:=\frac{1}{N}\sum_{i=1}^{N}\\{\ddot{u}(\hat{\bm{\theta}}^{T}_{MLE}\bm{x}_{i})\bm{x}_{i}\bm{x}^{T}_{i}[\dot{\psi}(u(\hat{\bm{\theta}}^{T}_{MLE}\bm{x}_{i}))-y_{i}]+\ddot{\psi}(u(\hat{\bm{\theta}}^{T}_{MLE}\bm{x}_{i}))\dot{u}^{2}(\hat{\bm{\theta}}^{T}_{MLE}\bm{x}_{i})\bm{x}_{i}\bm{x}^{T}_{i}\\}$

goes to a positive definite matrix in probability.

######  H 5.

Require that the full sample covariates have finite $6$-th order moments,
i.e., $E||\bm{x}_{1}||^{6}\leq\infty$.

######  H 6.

Assume $N^{-2}\sum_{i=1}^{N}||\bm{x}_{i}||^{s}/\phi_{i}=O_{P}(1)$ for $s=2,4$.

######  H 7.

For $\delta=0$ and some $\delta>0$, assume

$\frac{1}{N^{2+\delta}}\sum_{i=1}^{N}\frac{|y_{i}-\dot{\psi}_{i}(u(\hat{\bm{\theta}}^{T}_{MLE}\bm{x}_{i}))|^{2+\delta}||\dot{u}(\hat{\bm{\theta}}^{T}_{MLE}\bm{x}_{i})\bm{x}_{i}||^{2+\delta}}{\phi^{1+\delta}_{i}}=O_{P}(1).$

Assumptions H1 and H2 ensure that $\mbox{E}[y_{i}|\bm{x}_{i}]<\infty$ for all
$i$, and was used by [10] when investigating the impact of survey sampling
with unequal inclusion probabilities on (stochastic) gradient descent-based-
estimation methods in Big data problems. H2 defines an admissible set which is
required for ensuring consistency of estimates for model parameters for GLMs,
see [15]. Assumption H4 ensures (asymptotically) that $\bm{J_{X}}$, the
observed Fisher information, is defined for the given model, and that it is
positive definite and therefore non-singular. To obtain a Bahadur
representation, that is often useful in determining the asymptotic properties
of statistical estimators of the sub-sampled estimator, Assumptions H3 and H5
are needed. H6 and H7 are moment conditions on covariates and sub-sampling
probabilities. Assumption H7 is required by the Lindeberg-Feller central limit
theorem, see [35].

Based on the above assumptions, the following theorems were proved by [3]. The
first theorem proves that the estimator from the sub-sampling algorithm is
consistent where $||\tilde{\bm{\theta}}-\hat{\bm{\theta}}_{MLE}||$ can be made
small with a large sub-sample size $r$, and $\hat{\bm{\theta}}_{MLE}$ is the
maximum likelihood estimator of $\bm{\theta}$ based on $F_{N}$. The second
theorem proves that the approximation error,
$\tilde{\bm{\theta}}-\hat{\bm{\theta}}_{MLE}$, given $F_{N}$ is approximately
asymptotically Normally distributed with mean zero and variance $\bm{V}$.

######  Theorem 1.

If H1 to H7 hold then, as $N\rightarrow\infty$ and $r\rightarrow\infty$,
$\tilde{\bm{\theta}}$ is consistent to $\hat{\bm{\theta}}_{MLE}$ in
conditional probability given $F_{N}$. Moreover, the rate of convergence is
$r^{-\frac{1}{2}}$. That is, with probability approaching one, for any
$\epsilon>0$, there exist finite $\Delta_{\epsilon}$ and $r_{\epsilon}$ such
that

$P(||\tilde{\bm{\theta}}-\hat{\bm{\theta}}_{MLE}||\geq
r^{-\frac{1}{2}}\Delta_{\epsilon}|F_{N})<\epsilon,\forall r>r_{\epsilon}.$

######  Theorem 2.

If H1 to H7 hold, then as $N\rightarrow\infty$ and $r\rightarrow\infty$,
conditional on $F_{N}$ in probability,

$\bm{V}^{-1/2}(\tilde{\bm{\theta}}-\hat{\bm{\theta}}_{MLE})\rightarrow
N(\bm{0},\bm{I}),$

in distribution, where
$\bm{V}=\bm{J}^{-1}_{\bm{X}}\bm{V}_{c}\bm{J}^{-1}_{\bm{X}}=O_{P}(r^{-1})$ and

$\bm{V}_{c}=\frac{1}{rN^{2}}\sum_{i=1}^{N}\frac{\\{y_{i}-\dot{\psi}(u(\hat{\bm{\theta}}^{T}_{MLE}\bm{x}_{i}))\\}^{2}\dot{u}^{2}(\hat{\bm{\theta}}^{T}_{MLE}\bm{x}_{i})\bm{x}_{i}\bm{x}^{T}_{i}}{\phi_{i}}.$

When applying the general sub-sampling algorithm, it may not be clear how to
appropriately choose $\bm{\phi}$ depending upon the goal of the analysis
(e.g., parameter estimation, response prediction, etc). To address this, [3]
proposed determining $\bm{\phi}$ based on an optimality criterion from
experimental design (specifically $A$-optimality and $L$-optimality), and this
led to the proposal of Theorems 3 and 4 (below). Such theorems provide the
optimal choice for the sub-sampling probabilities to minimise the asymptotic
mean squared error of $\tilde{\bm{\theta}}$ (or $\mbox{tr}(\bm{V})$) and
$\bm{J_{X}}\tilde{\bm{\theta}}$ (or $\mbox{tr}(\bm{V}_{c})$), denoted as
$\bm{\phi}^{mMSE}$ and $\bm{\phi}^{mV_{c}}$, respectively. We note that such
sub-sampling probabilities are conditional on the assumption of a model being
appropriate to describe the Big data.

######  Theorem 3.

The sub-sampling strategy is $A$-optimal if the sub-sampling probability is
chosen such that

$\phi^{mMSE}_{i}=\frac{|y_{i}-\dot{\psi}(u(\hat{\bm{\theta}}^{T}_{MLE}\bm{x}_{i}))|\,||\bm{J}^{-1}_{\bm{X}}\dot{u}(\hat{\bm{\theta}}^{T}_{MLE}\bm{x}_{i})\bm{x}_{i}||}{\sum_{j=1}^{N}|y_{j}-\dot{\psi}(u(\hat{\bm{\theta}}^{T}_{MLE}\bm{x}_{j}))|\,||\bm{J}^{-1}_{\bm{X}}\dot{u}(\hat{\bm{\theta}}^{T}_{MLE}\bm{x}_{j})\bm{x}_{j}||},i=1,\ldots,N.$
(2)

Optimal sub-sampling probabilities from Theorem 3 can be computationally
expensive as they require the calculation of $\bm{J}^{-1}_{\bm{X}}$. Hence,
the linear- or $L$\- optimality criterion was proposed by [3], which minimises
the asymptotic mean squared error of $\bm{J_{X}}\tilde{\bm{\theta}}$ or
$\mbox{tr}(\bm{J_{X}}\bm{V}\bm{J_{X}})$ which corresponds to minimising the
variance of a linear combination of the parameters [5]. This led to the
following theorem.

######  Theorem 4.

The sub-sampling strategy is $L$-optimal if the sub-sampling probability is
chosen such that

$\phi^{mV_{c}}_{i}=\frac{|y_{i}-\dot{\psi}(u(\hat{\bm{\theta}}^{T}_{MLE}\bm{x}_{i}))|\,||\dot{u}(\hat{\bm{\theta}}^{T}_{MLE}\bm{x}_{i})\bm{x}_{i}||}{\sum_{j=1}^{N}|y_{j}-\dot{\psi}(u(\hat{\bm{\theta}}^{T}_{MLE}\bm{x}_{j}))|\,||\dot{u}(\hat{\bm{\theta}}^{T}_{MLE}\bm{x}_{j})\bm{x}_{j}||},i=1,\ldots,N.$
(3)

[40] applied Theorems 3 and 4 to obtain the following optimal sub-sampling
probabilities for logistic regression:

$\displaystyle\phi^{mMSE}_{i}=\frac{|y_{i}-\pi_{i}|\,||\bm{J}^{-1}_{\bm{X}}\bm{x}_{i}||}{\sum_{j=1}^{N}|y_{j}-\pi_{j}|\,||\bm{J}^{-1}_{\bm{X}}\bm{x}_{j}||},$
$\displaystyle\quad\phi^{mV_{c}}_{i}=\frac{|y_{i}-\pi_{i}|\,||\bm{x}_{i}||}{\sum_{j=1}^{N}|y_{j}-\pi_{j}|\,||\bm{x}_{j}||}$
(4)

with $y_{i}\in\\{0,1\\}$,
$\pi_{i}=\exp{(\hat{\bm{\theta}}^{T}_{MLE}\bm{x}_{i})}/(1+\exp{(\hat{\bm{\theta}}^{T}_{MLE}\bm{x}_{i})})$,
$\bm{J_{X}}=N^{-1}\sum_{h=1}^{N}\pi_{h}(1-\pi_{h})\bm{x}_{h}{\bm{x}_{h}}^{T}$
and $i=1,\ldots,N$.

[3] applied Theorems 3 and 4 to obtain the following optimal sub-sampling
probabilities for Poisson regression:

$\displaystyle\phi^{mMSE}_{i}=\frac{|y_{i}-\lambda_{i}|\,||\bm{J}^{-1}_{\bm{X}}\bm{x}_{i}||}{\sum_{j=1}^{N}|y_{j}-\lambda_{j}|\,||\bm{J}^{-1}_{\bm{X}}\bm{x}_{j}||},$
$\displaystyle\quad\phi^{mV_{c}}_{i}=\frac{|y_{i}-\lambda_{i}|\,||\bm{x}_{i}||}{\sum_{j=1}^{N}|y_{j}-\lambda_{j}|\,||\bm{x}_{j}||}$
(5)

with $y_{i}\in N_{0}$ or non-negative integers,
$\lambda_{i}=\exp{(\hat{\bm{\theta}}^{T}_{MLE}\bm{x}_{i})}$,
$\bm{J_{X}}=N^{-1}\sum_{h=1}^{N}\exp{(\hat{\bm{\theta}}^{T}_{MLE}\bm{x}_{h})}\bm{x}_{h}{\bm{x}_{h}}^{T}$
and $i=1,\ldots,N$.

Unfortunately, in practice, the optimal sub-sampling probabilities
$\bm{\phi}^{mMSE}$ and $\bm{\phi}^{mV_{c}}$ cannot be determined as they
depend on $\hat{\bm{\theta}}_{MLE}$. To address this, based on Theorems 1 and
2, [3] proposed a two stage sub-sampling strategy where an initial random
sample of the Big data is used to estimate $\hat{\bm{\theta}}_{MLE}$; an
estimate which is then used with the results from Theorems 3 and 4 to provide
estimates of optimal sub-sampling probabilities. Such an approach is thus
termed a two stage sub-sampling approach, and is outlined in the next section.

###  2.4 Optimal sub-sampling algorithm for GLMs

Theorems 1 to 4 provide a theoretical basis for the optimal sub-sampling
algorithm of [3]. In this section, we describe this approach for GLMs which is
outlined in Algorithm 2. This is a two stage algorithm where the general sub-
sampling algorithm is initially applied. For this, the sub-sampling
probabilities could be $1/N$ for all observations or may be based on a
stratified sampling technique. Based on this initial sub-sample, estimates
$\tilde{\bm{\theta}}$ are obtained. Then, optimal sub-sampling probabilities
are estimated using results from Theorems 3 and 4, where
$\hat{\bm{\theta}}_{MLE}$ is approximated by $\tilde{\bm{\theta}}$.

1

Stage 1

2 Random Sampling : Assign $\bm{\phi}=(\phi_{1},\ldots,\phi_{N})$ for $F_{N}$.
For example, in logistic regression $\phi_{i}=\phi^{prop}$ represents
proportional sub-sampling probabilities that are based on the response
composition, and $\phi_{i}=N^{-1}$ which is uniform sub-sampling probabilities
for Poisson regression.

According to $\bm{\phi}$ draw a random sub-sample of size $r_{0}$,
$S_{r_{0}}=\\{\bm{x}^{r_{0}}_{l},y^{r_{0}}_{l},\phi^{r_{0}}_{l}\\}_{l=1}^{r_{0}}=(\bm{X}^{r_{0}},\bm{y}^{r_{0}},\bm{\phi}^{r_{0}})$.

Estimation: Based on $S_{r_{0}}$, find:
$\displaystyle\tilde{\bm{\theta}}^{r_{0}}=\operatorname*{arg\,max}_{\bm{\theta}}~{}\log{L(\bm{\theta}|\bm{X}^{r_{0}},\bm{y}^{r_{0}},\bm{\phi}^{r_{0}})}\equiv\operatorname*{arg\,max}_{\bm{\theta}}~{}\frac{1}{r_{0}}\sum_{l=1}^{r_{0}}\Big{[}\frac{y^{r_{0}}_{l}u(\bm{\theta}^{T}\bm{x}^{r_{0}}_{l})-\psi(u(\bm{\theta}^{T}\bm{x}^{r_{0}}_{l}))}{\phi^{r_{0}}_{l}}\Big{]}.$

3 Stage 2

4 Optimal sub-sampling probability : Estimate $\bm{\phi}^{mMSE}$ or
$\bm{\phi}^{mV_{c}}$ from Theorems 3 and 4 using $\tilde{\bm{\theta}}^{r_{0}}$
in place of $\hat{\bm{\theta}}_{MLE}$.

5 Optimal sub-sampling and estimation: Based on the estimated sub-sampling
probabilities, draw a sub-sample $S_{r}$ completely at random (with
replacement) of size $r$ from $F_{N}$, such that
$S_{r}=\\{\bm{x}^{r}_{l},y^{r}_{l},\phi^{r}_{l}\\}_{l=1}^{r}=(\bm{X}^{r},\bm{y}^{r},\bm{\phi}^{r})$.

Form $S_{r_{0}+r}$ by combining $S_{r_{0}}$ and $S_{r}$, and obtain:
$\displaystyle\tilde{\bm{\theta}}=$
$\displaystyle\operatorname*{arg\,max}_{\bm{\theta}}~{}\log{L(\bm{\theta}|\bm{X}^{r_{0}},\bm{y}^{r_{0}},\bm{\phi}^{r_{0}},\bm{X}^{r},\bm{y}^{r},\bm{\phi}^{r})}$
$\displaystyle\equiv$
$\displaystyle\operatorname*{arg\,max}_{\bm{\theta}}~{}\frac{1}{r_{0}+r}\Bigg{[}\sum_{k\in\\{r_{0},r\\}}\sum_{l=1}^{k}\frac{y^{k}_{l}u(\bm{\theta}^{T}\bm{x}^{k}_{l})-\psi(u(\bm{\theta}^{T}\bm{x}^{k}_{l}))}{\phi^{k}_{l}}\Bigg{]}.$

Output : $\tilde{\bm{\theta}}$ and $S_{r_{0}+r}$.

Algorithm 2 Two stage optimal sub-sampling algorithm for GLMs.

The first stage of the algorithm is as shown in Algorithm 1. In the second
stage, $r\geq r_{0}$ data points are sampled with replacement from the Big
data through the optimal sub-sampling probabilities (estimated based on the
model parameters obtained from stage one). These two sub-sampled data sets are
then combined to yield a single data set, and this data set is used to
estimate the model parameters. Once these estimates have been obtained,
$\tilde{\bm{V}}$, the variance-covariance matrix of $\tilde{\bm{\theta}}$, can
be estimated via
$\tilde{\bm{V}}=\tilde{\bm{J}}^{-1}_{\bm{X}}\tilde{\bm{V}}_{c}\tilde{\bm{J}}^{-1}_{\bm{X}}$,
where

$\displaystyle\tilde{\bm{J_{X}}}=$
$\displaystyle\frac{1}{N(r_{0}+r)}\Bigg{[}\sum_{k\in\\{r_{0},r\\}}\sum_{l=1}^{k}\frac{\ddot{u}(\tilde{\bm{\theta}}^{T}\bm{x}^{k}_{l})\bm{x}^{k}_{l}(\bm{x}^{k}_{l})^{T}[\ddot{\psi}(u(\tilde{\bm{\theta}}^{T}\bm{x}^{k}_{l}))-y^{k}_{l}]+\ddot{\psi}(u(\tilde{\bm{\theta}}^{T}\bm{x}^{k}_{l}))\dot{u}^{2}(\tilde{\bm{\theta}}^{T}\bm{x}^{k}_{l})\bm{x}^{k}_{l}(\bm{x}^{k}_{l})^{T}}{\phi^{k}_{l}}\Bigg{]}$
and
$\displaystyle\quad\tilde{\bm{V}_{c}}=\frac{1}{N^{2}(r_{0}+r)^{2}}\Bigg{[}\sum_{k\in\\{r_{0},r\\}}\sum_{l=1}^{k}\frac{[y^{k}_{l}-\dot{\psi}(u(\tilde{\bm{\theta}}^{T}\bm{x}^{k}_{l}))]^{2}\dot{u}^{2}(\tilde{\bm{\theta}}^{T}\bm{x}^{k}_{l})\bm{x}^{k}_{l}(\bm{x}^{k}_{l})^{T}}{(\phi^{k}_{l})^{2}}\Bigg{]}.$
(6)

Specific versions of the above algorithm under the logistic and Poisson
regression models are given in Appendix B in Algorithms 4 and 5, respectively.

###  2.5 Limitations of the optimal sub-sampling algorithm

The above optimal sub-sampling approach has a number of limitations. One such
limitation is the computational expense involved in obtaining the optimal sub-
sampling probabilities, as these need to be found for each data point in the
Big data set. To address this, [21] introduced a faster two stage sub-sampling
procedure for GLMs using the Johnson-Lindenstrauss Transform (JLT) and sub-
sampled Randomised Hadamard Transform (SRHT), which are techniques to downsize
matrix volume. Another limitation is that, as the approximated optimal sub-
sampling probabilities are proportional to
$|y_{i}-\dot{\psi}(u(\tilde{\bm{\theta}}^{T}\bm{x}_{i}))|$, an observation
with $y_{i}\approx\dot{\psi}(u(\tilde{\bm{\theta}}^{T}\bm{x}_{i}))$ has a near
zero probability of being selected and data-points with
$y_{i}=\dot{\psi}(u(\tilde{\bm{\theta}}^{T}\bm{x}_{i}))$ will be never sub-
sampled. [3] proposed to resolve this by introducing $\epsilon$ (a small
positive value, e.g., $10^{-6}$) to constrain the optimal sub-sampling
probabilities by replacing
$[y_{i}-\dot{\psi}(u(\tilde{\bm{\theta}}^{T}\bm{x}_{i}))]$ with
$\max{(|y_{i}-\dot{\psi}(u(\tilde{\bm{\theta}}^{T}\bm{x}_{i}))|,\epsilon)}$
which ensures such data points have a non-zero probability of being selected.
Lastly, one of the major limitations of the approach that has not been
addressed previously is the inherent assumption that the Big data can be
appropriately described by a given model. That is, the sub-sampling
probabilities are evaluated based on an assumed model, and they are only
optimal for this model. We suggest that this is a substantial limitation as
specifying such a model in practice can be difficult. This motivates the
development of methods that yield sub-sampling probabilities that are robust
to the choice of model, and our proposed approach for this is outlined next.

##  3 Model robust optimal sub-sampling method

In order to apply the above two stage sub-sampling approach, optimal sub-
sampling probabilities need to be evaluated, and these are based on a model
that is assumed to appropriately describe the data. In practice, determining
such a model may be difficult, and there could be a variety of models that
appropriately describe the data. Hence, a sub-sampling approach that provides
robustness to the choice of model is desirable. For this, we propose to
consider a set of $Q$ models which can be constructed to encapsulate a variety
of scenarios that may be observed within the data. For each model in this set,
define model probabilities $\alpha_{q}$ for $q=1,\ldots,Q$ such that
$\sum_{q=1}^{Q}\alpha_{q}=1$, which represents our a priori belief about the
appropriateness of each model. Denote the model matrix for the $q$-th model as
$\bm{X}_{q}=h_{q}(\bm{X}_{0})$, i.e., some function of the data matrix
$\bm{X}_{0}$. To apply a sub-sampling approach for this model set, sub-
sampling probabilities are needed, and they should be constructed such that
the resulting data is expected to address the analysis aim, regardless of
which model is actually preferred for the Big data. For this purpose, we
propose to form these sub-sampling probabilities by taking a weighted average
(based on $\alpha_{q})$ of the sub-sampling probabilities that would be
obtained for each model (singularly). This is the basic approach of our model
robust optimal sub-sampling algorithm, for which further details are provided
below, including a theoretical basis for constructing the sub-sampling
probabilities in this way.

###  3.1 Properties of model robust optimal sub-sampling algorithm

Updating the notation of the model matrix from $\bm{X}$ to $\bm{X}_{q}$ for
the $q$-th model subsequently leads to analogous definitions for
$\bm{x}_{qi}$, $\bm{\theta}_{q}$, $\hat{\bm{\theta}}_{qMLE}$,
$\tilde{\bm{\theta}}_{q}$, ${\bm{J_{X}}}_{q}$, $\bm{V}_{q}$ and
${\bm{V}_{q}}_{c}$. Assuming H1 to H7 hold for each of the $q$ models,
Theorems 1 and 2 apply straightforwardly to each of the models. Extending the
ideas from Section 2.3, optimal sub-sampling probabilities can be selected
based on certain optimality criteria to, for example, ensure efficient
estimates of parameters across the $Q$ models. This leads to the following
theorems.

######  Theorem 5.

For a set of $Q$ models with model probability $\alpha_{q}$ for the $q$-th
model, $q=1,\ldots Q$, if the sub-sampling probabilities are selected as
follows:

$\phi^{mMSE}_{i}=\sum_{q=1}^{Q}\alpha_{q}\frac{|y_{i}-\dot{\psi}(u(\hat{\bm{\theta}_{q}}_{MLE}^{T}\bm{x}_{qi}))|\,||\bm{J}^{-1}_{\bm{X}_{q}}\dot{u}(\hat{\bm{\theta}_{q}}_{MLE}^{T}\bm{x}_{qi})\bm{x}_{qi}||}{\sum_{j=1}^{N}|y_{j}-\dot{\psi}(u(\hat{\bm{\theta}_{q}}_{MLE}^{T}\bm{x}_{qj}))|\,||\bm{J}^{-1}_{\bm{X}_{q}}\dot{u}(\hat{\bm{\theta}_{q}}_{MLE}^{T}\bm{x}_{qj})\bm{x}_{qj}||},$
(7)

$i=1,\ldots,N$, and $\sum_{q=1}^{Q}\alpha_{q}=1$, then
$\sum_{q=1}^{Q}\alpha_{q}\mbox{tr}(\bm{V}_{q})$ attains its minimum.

The proof of Theorem 5 is available in Appendix A, which is an extension of
the proof of Theorem 3 in [3].

######  Theorem 6.

For a set of $Q$ models with model probability $\alpha_{q}$ for the $q$-th
model, $q=1,\ldots Q$, if the sub-sampling probabilities are selected as
follows:

${\phi}_{i}^{mV_{c}}=\sum_{q=1}^{Q}\alpha_{q}\frac{|y_{i}-\dot{\psi}(u(\hat{\bm{\theta}_{q}}_{MLE}^{T}\bm{x}_{qi}))|\,||\dot{u}(\hat{\bm{\theta}_{q}}_{MLE}^{T}\bm{x}_{qi})\bm{x}_{qi}||}{\sum_{j=1}^{N}|y_{j}-\dot{\psi}(u(\hat{\bm{\theta}_{q}}_{MLE}^{T}\bm{x}_{qj}))|\,||\dot{u}(\hat{\bm{\theta}_{q}}_{MLE}^{T}\bm{x}_{qj})\bm{x}_{qj}||},$
(8)

$i=1,\ldots,N$, and $\sum_{q=1}^{Q}\alpha_{q}=1$, then
$\sum_{q=1}^{Q}\alpha_{q}\mbox{tr}({\bm{V}_{q}}_{c})$ attains its minimum.

The proof of Theorem 6 follows straightforwardly from the proof of Theorem 5.

The above theorems can be applied to obtain optimal sub-sampling probabilities
under the logistic regression model as follows:

$\displaystyle{\phi}^{mMSE}_{i}=\sum_{q=1}^{Q}\alpha_{q}\frac{|y_{i}-{\pi_{q}}_{i}|\,||\bm{J}^{-1}_{\bm{X}_{q}}\bm{x}_{qi}||}{\sum_{j=1}^{N}|y_{j}-{\pi_{q}}_{j}|\,||\bm{J}^{-1}_{\bm{X}_{q}}\bm{x}_{qj}||},$
$\displaystyle\quad{\phi}^{mV_{c}}_{i}=\sum_{q=1}^{Q}\alpha_{q}\frac{|y_{i}-{\pi_{q}}_{i}|\,||\bm{x}_{qi}||}{\sum_{j=1}^{N}|y_{j}-{\pi_{q}}_{j}|\,||\bm{x}_{qj}||}$
(9)

with $y_{i}\in\\{0,1\\}$,
${\pi_{q}}_{i}=\exp{(\hat{\bm{\theta}}_{qMLE}^{T}\bm{x}_{qi})}/(1+\exp{(\hat{\bm{\theta}}_{qMLE}^{T}\bm{x}_{qi})}),$
$\bm{J}_{\bm{X}_{q}}=N^{-1}\sum_{h=1}^{N}{\pi_{q}}_{h}(1-{\pi_{q}}_{h})\bm{x}_{qh}(\bm{x}_{qh})^{T}$
and $i=1,\ldots,N$.

Similarly, optimal sub-sampling probabilities under the Poisson regression
model can be obtained as follows:

$\displaystyle{\phi}^{mMSE}_{i}=\sum_{q=1}^{Q}\alpha_{q}\frac{|y_{i}-{\lambda_{q}}_{i}|\,||\bm{J}^{-1}_{\bm{X}_{q}}\bm{x}_{qi}||}{\sum_{j=1}^{N}|y_{j}-{\lambda_{q}}_{j}|\,||\bm{J}^{-1}_{\bm{X}_{q}}\bm{x}_{qj}||},$
$\displaystyle\quad\phi^{mV_{c}}_{i}=\sum_{q=1}^{Q}\alpha_{q}\frac{|y_{i}-{\lambda_{q}}_{i}|\,||\bm{x}_{qi}||}{\sum_{j=1}^{N}|y_{j}-{\lambda_{q}}_{j}|\,||\bm{x}_{qj}||}$
(10)

with $y_{i}\in N_{0}$ or non-negative integers,
${\lambda_{q}}_{i}=\exp{(\hat{\bm{\theta}}_{qMLE}^{T}\bm{x}_{qi})},$
$\bm{J}_{\bm{X}_{q}}=N^{-1}\sum_{h=1}^{N}{\lambda_{q}}_{h}\bm{x}_{qh}(\bm{x}_{qh})^{T}$
and $i=1,\ldots,N$.

As noted in Section 2.3 and by [3], these model robust optimal sub-sampling
probabilities are based on the maximum likelihood estimator found by
considering the whole Big data set. Hence, a two stage procedure similar to
Algorithm 2 is proposed for model robust sub-sampling, and this is outlined
next.

###  3.2 Model robust optimal sub-sampling algorithm for GLMs

The two stage model robust optimal sub-sampling algorithm for GLMs is
presented in Algorithm 3, where the estimates of sub-sampling probabilities
based on the results of Theorems 5 and 6 are used. The specific algorithms for
logistic and Poisson regression are available in Appendix B as Algorithms 6
and 7, respectively.

1

Stage 1

2 Random Sampling : Assign $\bm{\phi}=(\phi_{1},\ldots,\phi_{N})$ for $F_{N}$.
For example, in logistic regression $\phi_{i}=\phi^{prop}$ represents
proportional sub-sampling probabilities that are based on the response
composition, and $\phi_{i}=N^{-1}$ which is uniform sub-sampling probabilities
for Poisson regression.

According to $\bm{\phi}$ draw random sub-samples of size $r_{0}$, such that
$S_{r_{0}}=\\{h_{q}({\bm{x}_{0}}^{r_{0}}_{l}),y^{r_{0}}_{l},\phi^{r_{0}}_{l}\\}_{l=1}^{r_{0}}=(h_{q}(\bm{X}_{0}^{r_{0}}),\bm{y}^{r_{0}},\bm{\phi}^{r_{0}})$
for $q=1,\ldots,Q$.

Estimation: For $q=1,\ldots,Q$ and $S_{r_{0}}$, find:
$\displaystyle\tilde{\bm{\theta}^{r_{0}}_{q}}=\operatorname*{arg\,max}_{\bm{\theta}_{q}}~{}\log{L(\bm{\theta}_{q}|\bm{X}^{r_{0}}_{q},\bm{y}^{r_{0}},\bm{\phi}^{r_{0}})}\equiv\operatorname*{arg\,max}_{\bm{\theta}_{q}}~{}\frac{1}{r_{0}}\sum_{l=1}^{r_{0}}\Big{[}\frac{y^{r_{0}}_{l}u(\bm{\theta}_{q}^{T}\bm{x}^{r_{0}}_{ql})-\psi(u(\bm{\theta}_{q}^{T}\bm{x}^{r_{0}}_{ql})}{\phi^{r_{0}}_{l}}\Big{]}.$

3 Stage 2

4 Optimal sub-sampling probability: Estimate optimal sub-sampling
probabilities $\bm{\phi}^{mMSE}$ or $\bm{\phi}^{mV_{c}}$ using the results
from Theorem 5 or 6, with $\hat{\bm{\theta}}_{qMLE}$ replaced by
$\tilde{\bm{\theta}}^{r_{0}}_{q}$.

5 Optimal sub-sampling and estimation: Based on $\bm{\phi}^{mMSE}$ or
$\bm{\phi}^{mV_{c}}$, draw sub-sample of size $r$ from $F_{N}$, such that
$S_{r}=\\{h_{q}({\bm{x}_{0}}^{r}_{l}),\bm{y}^{r}_{l},\bm{\phi}^{r}_{l}\\}_{l=1}^{r}=(h_{q}(\bm{X}_{0}^{r}),\bm{y}^{r},\bm{\phi}^{r})$
for $q=1,\ldots,Q$.

Combine $S_{r_{0}}$ and $S_{r}$ to form $S_{(r_{0}+r)}$, and obtain :
$\displaystyle\tilde{\bm{\theta}_{q}}=$
$\displaystyle\operatorname*{arg\,max}_{\bm{\theta}_{q}}~{}\log{\big{(}L(\bm{\theta}_{q}|\bm{X}^{r_{0}}_{q},\bm{y}^{r_{0}},\bm{\phi}^{r_{0}},\bm{X}^{r}_{q},\bm{y}^{r},\bm{\phi}^{r})\big{)}}$
$\displaystyle\equiv$
$\displaystyle\operatorname*{arg\,max}_{\bm{\theta}_{q}}~{}\frac{1}{r_{0}+r}\Bigg{[}\sum_{k\in\\{r_{0},r\\}}\sum_{l=1}^{k}\frac{y^{k}_{l}u(\bm{\theta}^{T}_{q}\bm{x}^{k}_{ql})-\psi(u(\bm{\theta}^{T}_{q}\bm{x}^{k}_{ql}))}{\phi^{k}_{l}}\Bigg{]}.$

Output: $\tilde{\bm{\theta}}_{q}$ and $S_{(r_{0}+r)}$ for $q=1,\ldots,Q$.

Algorithm 3 Two stage model robust optimal sub-sampling algorithm for GLMs.

Similar to Algorithm 2, the first stage entails randomly sub-sampling $F_{N}$
(with replacement), and estimating model parameters for each of the $Q$
models. Based on these estimated parameters, model specific sub-sampling
probabilities are obtained, and these are combined based on $\alpha_{q}$,
forming model robust optimal sub-sampling probabilities. Subsequently, $r\geq
r_{0}$ data points are sampled from $F_{N}$. The two sub-samples are then
combined, and each of the $Q$ models are fitted (separately) based on the
weighted log-likelihood function, which should yield efficient estimates of
parameters across each of the models. Similar to estimating the variance-
covariance matrix of a single model (Equation (2.4)), $\tilde{\bm{V}}_{q}$,
the variance-covariance matrix of $\tilde{\bm{\theta}_{q}}$ can be estimated,
for $q=1,\ldots,Q$.

In the following section, the proposed model robust optimal sub-sampling
method and the current optimal sub-sampling method are assessed through a
simulation study and real world scenarios.

##  4 Applications of optimal sub-sampling algorithms

In this section, a simulation study and two real world applications are used
to assess the performance of the proposed model robust optimal sub-sampling
algorithm (Algorithm 3) compared to the optimal sub-sampling algorithm
(Algorithm 2), and random sampling. The main results are presented in this
section with some results presented in the Appendix. The simulation study and
real world applications were coded in the R statistical programming language
with the help of Rstudio IDE, and our scripts are available through Github.
Appendix C provides specific Github hyperlinks to the code repositories for
the simulation study and real world applications.

###  4.1 Simulation study design

To explore the performance of our model robust sub-sampling approach, a
simulation study was constructed based on the logistic and Poisson regression
models. For each case, a set of $Q=4$ models were assumed based on [34], and
this set is summarised in Table 1. For each model, $F_{N}$ was constructed, by
assuming a distribution for the covariates and corresponding response. The
performance of the three sampling methods were then compared for each $F_{N}$
through evaluating six scenarios: 1) random sampling to estimate the
parameters of the data generating model; 2) optimal sub-sampling under the
data generating model – this simulates the case where an appropriate model was
assumed for describing the Big data; 3)-5) optimal sub-sampling under
alternative models i.e., optimal sub-sampling to estimate the parameters of
the data generating model but where samples were obtained based on an
alternative model – this simulates undertaking optimal sub-sampling based on a
‘wrong’ model; 6) model robust optimal sub-sampling, with $\alpha_{q}=1/4$
used to estimate the parameters of the data generating model, for
$q=1,\ldots,4$. To quantitatively compare each approach, the simulated Mean
Squared Error (SMSE) was used. This was evaluated as follows:

$SMSE(\bm{\theta})=\frac{1}{M}{\sum_{m=1}^{M}\sum_{n=1}^{p+q}(\tilde{\theta}_{nm}-\theta_{n})^{2}},$
(11)

where $M$ is the number of simulations, $p+q$ is the number of parameters,
$\theta_{n}$ is the $n$-th underlying parameter of the data generating model
and $\tilde{\theta}_{nm}$ is the estimate of this parameter from the $m$-th
simulation.

In addition, the determinant of the observed Fisher information matrix (i.e.,
the inverse of the variance-covariance matrix $\bm{V}$) was also compared
across sub-sampling approaches, with the largest of such values being
preferred. For each simulation, $N=10000$, $r_{0}=100$,
$r=100,200,\ldots,1400$ and $M=1000$.

Table 1: Model set assumed for simulation study. Model set
---
$\theta_{0}+\theta_{1}\bm{x}_{1}+\theta_{2}\bm{x}_{2}$
$\theta_{0}+\theta_{1}\bm{x}_{1}+\theta_{2}\bm{x}_{2}+\theta_{3}\bm{x}^{2}_{1}$
$\theta_{0}+\theta_{1}\bm{x}_{1}+\theta_{2}\bm{x}_{2}+\theta_{3}\bm{x}^{2}_{2}$
$\theta_{0}+\theta_{1}\bm{x}_{1}+\theta_{2}\bm{x}_{2}+\theta_{3}\bm{x}^{2}_{1}+\theta_{4}\bm{x}^{2}_{2}$

Additionally, we only consider the $A$-optimal ($mMSE$) sub-sampling strategy
in these simulations, as [40, 3] indicated that this approach generally
outperformed the $L$-optimal ($mV_{c}$) strategy. We explore both optimality
criteria in the real-world applications.

####  4.1.1 Logistic regression

Following [40], covariate data for the logistic regression model were
simulated based on two distributions: Exponential ($\lambda$) and Multivariate
Normal ($\bm{\mu},\bm{\Sigma}$). The values of
($\lambda,\bm{\mu},\bm{\Sigma}$) and $\bm{\theta}$ are given in Table 2 for
each data generating model. For all models, the first element of $\bm{\theta}$
is the value of the parameter for the intercept. While $\lambda$, $\bm{\mu}$
and $\bm{\Sigma}$ were selected arbitrarily, $\bm{\theta}$ was selected such
that for $N=10000$, the data generating model was preferred over the
alternative models based on the Akaike Information Criterion [4].

Table 2: $\lambda,\bm{\mu},\bm{\Sigma}$ and $\bm{\theta}$ values are given to generate $x_{1},x_{2}$ through Exponential and Multivariate Normal distributions and form $F_{N}$ for each data generating logistic regression model. Covariates | Distributions
---|---
Exponential ($\lambda=\sqrt{3}$) |  Normal $\Big{(}\bm{\mu}=[0,0]$,$\bm{\Sigma}=\begin{bmatrix}1.5&0\\\ 0&1.5\end{bmatrix}\Big{)}$
$\bm{x}_{1},\bm{x}_{2}$ | $\bm{\theta}=[-2,\hphantom{-}1.5,\hphantom{-}0.3]$ | $\bm{\theta}=[-1,\hphantom{-}0.5,\hphantom{-}0.1]$
$\bm{x}_{1},\bm{x}_{2},\bm{x}^{2}_{1}$ | $\bm{\theta}=[-2,\hphantom{-}1.7,-1.2,\hphantom{-}0.2]$ | $\bm{\theta}=[-1,\hphantom{-}0.5,\hphantom{-}0.5,\hphantom{-}0.7]$
$\bm{x}_{1},\bm{x}_{2},\bm{x}^{2}_{2}$ | $\bm{\theta}=[-2,-1.3,\hphantom{-}1.9,\hphantom{-}0.9]$ | $\bm{\theta}=[-1,\hphantom{-}0.5,\hphantom{-}0.5,\hphantom{-}0.3]$
$\bm{x}_{1},\bm{x}_{2},\bm{x}^{2}_{1},\bm{x}^{2}_{2}$ | $\bm{\theta}=[-2,\hphantom{-}1.9,\hphantom{-}1.9,\hphantom{-}0.9,\hphantom{-}0.7]$ | $\bm{\theta}=[-1,\hphantom{-}0.5,\hphantom{-}0.5,\hphantom{-}0.5,\hphantom{-}0.5]$
Figure 1: (a) Model information (larger is better) and (b) SMSE (smaller is
better) for the sub-sampling methods for the logistic regression model under
$mMSE$. Covariate data were generated from the Multivariate Normal
distribution.

Figure 1 provides summaries of the SMSE and average model information over all
models when the covariate data is generated through a Multivariate Normal
distribution. The SMSE and average model information indicate that, under
optimal sub-sampling for $mMSE$, the data generating model is typically
preferred within the model set. This is expected as it is the case where the
appropriate data generating model was correctly assumed to describe the Big
data. Of note, the proposed model robust approach performs similarly to the
optimal sub-sampling approach. Notable increases in the SMSE and decreases in
the model information are observed when the incorrect model is considered for
optimal sub-sampling compared to the data generating model. Overall, random
sampling tends to have the worst performance. Similar results were obtained
when the covariate data were generated through an Exponential distribution
(Appendix D).

####  4.1.2 Poisson regression

The simulation study based on Poisson regression was constructed similar to
the logistic regression case. In terms of generating covariate values, Uniform
and Multivariate Normal ($\bm{\mu},\bm{\Sigma}$) distributions were used, see
Table 3 and [3]. Values for $\bm{\theta}$ were selected as described above,
and are given in this table.

Table 3: $\bm{\mu},\bm{\Sigma}$ and $\bm{\theta}$ values are given to generate $x_{1},x_{2}$ through uniform and Multivariate Normal distributions and form $F_{N}$ for each data generating Poisson regression model. Covariates | Distributions
---|---
Uniform |  Normal $\Big{(}\bm{\mu}=[0,0]$,$\bm{\Sigma}=\begin{bmatrix}1&0\\\ 0&1\end{bmatrix}\Big{)}$
$\bm{x}_{1},\bm{x}_{2}$ | $\bm{\theta}=[1,\hphantom{-}0.3,\hphantom{-}0.3]$ | $\bm{\theta}=[1,\hphantom{-}0.5,\hphantom{-}0.1]$
$\bm{x}_{1},\bm{x}_{2},\bm{x}^{2}_{1}$ | $\bm{\theta}=[1,\hphantom{-}0.7,-0.5,\hphantom{-}0.4]$ | $\bm{\theta}=[1,-0.3,\hphantom{-}0.3,\hphantom{-}0.1]$
$\bm{x}_{1},\bm{x}_{2},\bm{x}^{2}_{2}$ | $\bm{\theta}=[1,-0.4,\hphantom{-}0.5,\hphantom{-}0.3]$ | $\bm{\theta}=[1,\hphantom{-}0.3,-0.3,\hphantom{-}0.1]$
$\bm{x}_{1},\bm{x}_{2},\bm{x}^{2}_{1},\bm{x}^{2}_{2}$ | $\bm{\theta}=[1,\hphantom{-}0.5,\hphantom{-}0.5,-0.3,\hphantom{-}0.5]$ | $\bm{\theta}=[1,\hphantom{-}0.3,\hphantom{-}0.3,-0.3,\hphantom{-}0.1]$

SMSE and average model information over all models when the covariate data was
generated from a Multivariate Normal distribution and the response from
generated form a Poisson regression model are shown in Figure 2. Generally,
random sampling performs worst, while the proposed model robust approach and
the optimal sub-sampling method based on the data generating model perform the
best, and have similar SMSE and average model information values. Again, the
use of the optimal sub-sampling algorithm can lead to notable increases in
SMSE when the assumed model is incorrect.

Similar results were obtained when the covariate data was generated from a
uniform distribution (see Appendix D).

Figure 2: (a) Model information (larger is better) and (b) SMSE (smaller is
better) for the sub-sampling methods for the Poisson regression model under
$mMSE$. Covariate data were generated from the standardised Multivariate
Normal distribution.

###  4.2 Real world applications

The three sub-sampling methods are applied to analyse the “Skin segmentation"
and “New York City taxi fare" data under logistic and Poisson regression,
respectively. In the simulation study, the parameters were specified for the
data generating model. However, in real world applications these are unknown.
In this situation the sub-sampling methods cannot be compared as in Section
4.1. Instead, for every model, the simulated mean squared error (SMSE) under
each sub-sampling method is evaluated for various $r$ sub-sample sizes and $M$
simulations. This can be evaluated as follows:

$SSMSE(\hat{\bm{\theta}}_{MLE})=\sum_{q=1}^{Q}SMSE_{q}(\hat{\bm{\theta}}_{MLE}),$
(12)

where $SMSE_{q}(\hat{\bm{\theta}}_{MLE})$ is the simulated mean squared error
in Equation (11) with $\bm{\theta}$ replaced by $\hat{\bm{\theta}}_{MLE}$.

In the following real world examples, the set of $Q$ models includes the main
effects model, with intercept, and all possible combinations of quadratic
terms for continuous covariates. Again, these were constructed based on the
work of [34].

####  4.2.1 Identifying skin from colours in images

[30] considered the problem of identifying skin-like regions in images as part
of the complex task of performing facial recognition. For this purpose, [30]
collated the “Skin segmentation" data set, which consists of RGB (R-red,
G-green, B-blue) values of randomly sampled pixels from $N=245,057$ face
images (out of which $50,859$ are skin samples and $194,198$ are non-skin
samples) from various age groups, race groups and genders. [6, 7] applied
multiple supervised machine learning algorithms to classify if images are skin
or not based on the RGB colour data. In addition, [1] conducted the same
classification task for two different colour spaces, HSV (H-hue, S-saturation,
V-value) and YCbCr (Y-luma component, Cb-blue difference chroma component, Cr-
red difference chroma component), by transforming the RGB colour space.

We consider the same classification problem but use a logistic regression
model. Skin presence is denoted as one and skin absence is denoted as zero.
Each colour vector is scaled to have a mean of zero and a variance of one
(initial range was between $0-255$). To compare sub-sampling methods, we set
$r_{0}=200,r=200,300,\ldots,1800$ for the sub-samples, and construct a set of
$Q$ models by considering an intercept with main effects model with all
covariates (scaled colors red,green and blue) as the base model, and form all
alternative models by including different combinations of quadratic terms of
all covariates. This leads to $Q=8$. Each of these models is considered
equally likely a priori.

Figure 3 shows the SSMSE values over the $Q$ models for various sub-sample
sizes obtained by applying logistic regression to the “Skin segmentation"
data. The proposed model robust approach performs similarly to the optimal
sub-sampling method for $r=200$ and $r=300$. However, when the sample size
increases the optimal sub-sampling method performs poorly, and after $r=1300$
random sampling actually has lower SSMSE values than optimal sub-sampling. The
same is not true for our proposed model robust approach which has the lowest
values of SSMSE throughout the selected sample sizes under $mMSE$ and $mV_{c}$
optimal sub-sampling criteria. It is interesting that the optimal sub-sampling
approach based on $mMSE$ performed worse than random sampling in comparison to
the simulation study results. Upon investigating this, it was found that it
could potentially be explained by one of the models being particularly poor
(subsequent higher SMSE with increasing $r$ values) for describing the data,
and therefore led to inflated SSMSE values. Despite this, we note that the
model robust approach appears to perform well in general.

Figure 3: Summed SMSE over the available models for logistic regression
applied on the “Skin segmentation" data.

####  4.2.2 New York City taxi cab usage

New York City (NYC) taxi trip and fare information from 2009 onward,
consisting of over $170$ million records each year, are publicly available
courtesy of the New York City Taxi and Limousine Commission. Some analyses of
interest of these NYC taxi data include: a taxi driver’s decision process to
pick up a fair or cruise for customers which was modelled via a logistic
regression model [46]; taxi demand and how it is impacted by location, time,
demographic information, socioeconomic status and employment status which was
modelled via a multiple linear regression model [42]; and the dependence of
taxi supply and demand on location and time considered via the Poisson
regression model [43].

In our application, we are interested in how taxi usage varies with day of the
week (weekday/weekend), season (winter/spring/summer/autumn), fare amount, and
cash payment. The data used in our application is the “New York City taxi
fare” data for the year 2013, hosted by the University of Illinois Urbana
Champaign [13].

Each data point includes the number of rides recorded against the medallion
number (a license number provided for a motor vehicle to operate as a taxi)
($y$), weekday or not ($\bm{x_{1}}$), winter or not ($\bm{x_{2}}$), spring or
not ($\bm{x_{3}}$), summer or not ($\bm{x_{4}}$), summed fare amount in
dollars ($\bm{x_{5}}$), and the ratio of cash payment trips in terms of all
trips ($\bm{x_{6}}$). The continuous covariate $\bm{x_{5}}$ was scaled to have
a range of zero and one. Poisson regression was used to model the relationship
between the number of rides per medallion and these covariates.

Figure 4: Summed SMSE over the available models for Poisson regression applied
on the 2013 NYC taxi usage data.

For our study, three sub-sampling methods are compared for the analysis of the
taxi fare data, assigning $r_{0}=100$ and $r=200,\ldots,1900$ for sub-samples.
The model set consists of the main effects
model($\bm{x_{1}},\bm{x_{2}},\bm{x_{3}},\bm{x_{4}},\bm{x_{5}},\bm{x_{6}}$) and
all possible combinations of the quadratic terms of the continuous covariates
($\bm{x^{2}_{5}},\bm{x^{2}_{6}}$) which leads to $Q=4$. Each of these models
were considered equally likely a priori.

The SSMSE over the four models is shown in Figure 4. Our proposed model robust
approach outperforms the optimal sub-sampling method for almost all sample
sizes for both $mMSE$ and $mV_{c}$ strategy. Under $mV_{c}$, the model robust
approach actually initially performs worse than the optimal sub-sampling
approach but this is quickly reserved as $r$ is increased. Random sampling
performs the worst, suggesting that there is benefit in using targeted
sampling approaches (as proposed here) over random selection.

##  5 Discussion

In this article, we proposed a model robust optimal sub-sampling approach for
GLMs. This new approach extends the current optimal sub-sampling approach by
considering a set of models (rather than a single model) when determining the
sub-sampling probabilities. The exact formulation of these probabilities is
derived using Theorems 5 and 6. The robustness properties of this proposed
approach were demonstrated in a simulation study, and in two real-world
analysis problems, where it was shown to outperform optimal sub-sampling and
random sampling. Accordingly, we suggest that such an approach be considered
in future Big data analysis problems.

The main limitation of the proposed approach that could be addressed in future
research is extending the specification of the model set to a flexible class
of models. This could be, for example, through a generalised additive model
[17] or the inclusion of a discrepancy term in the linear predictor [20].
Another avenue of interest that could also be explored is reducing the model
set after stage one of the two-stage algorithm where models that clearly do
not appear to be appropriate for the data could be dropped. Both of these
extensions are planned for future research.

## References

  * [1] Ayad R Abbas and Ayat O Farooq  “Skin Detection using Improved ID3 Algorithm”  In _Iraqi Journal of Science_ 60.2, 2019, pp. 402–410  URL: https://ijs.uobaghdad.edu.iq/index.php/eijs/article/view/658
  * [2] Mingyao Ai, Fei Wang, Jun Yu and Huiming Zhang  “Optimal subsampling for large-scale quantile regression”  In _Journal of Complexity_ 62 Elsevier, 2021, pp. 101512  DOI: https://doi.org/10.1016/j.jco.2020.101512
  * [3] Mingyao Ai, Jun Yu, Huiming Zhang and HaiYing Wang  “Optimal subsampling algorithms for Big data regressions”  In _Statistica Sinica_ 31, 2021, pp. 749–772  DOI: https://doi.org/10.5705/ss.202018.0439
  * [4] Hirotugu Akaike  “A new look at the statistical model identification”  In _IEEE Transactions on Automatic Control_ 19.6 IEEE, 1974, pp. 716–723  DOI: https://doi.org/10.1109/TAC.1974.1100705
  * [5] Anthony Atkinson, Alexander Donev and Randall Tobias  “Optimum Experimental Designs, with SAS”, Oxford Statistical Science Series  OUP Oxford, 2007
  * [6] Rajen B Bhatt, Gaurav Sharma, Abhinav Dhall and Santanu Chaudhury  “Efficient Skin Region Segmentation Using Low Complexity Fuzzy Decision Tree Model”  In _2009 Annual IEEE India Conference_ , 2009, pp. 1–4  IEEE  DOI: https://doi.org/10.1109/INDCON.2009.5409447
  * [7] Bartosz Binias, Mariusz Frąckiewicz, Krzysztof Jaskot and Henryk Palus  “Pixel Classification for Skin Detection in Color Images”  In _Advanced Technologies in Practical Applications for National Security_ Springer International Publishing, 2018, pp. 87–99  DOI: https://doi.org/10.1007/978-3-319-64674-9_6
  * [8] Xiangyu Chang, Shao-Bo Lin and Yao Wang  “Divide and conquer local average regression”  In _Electronic Journal of Statistics_ 11.1 The Institute of Mathematical Statisticsthe Bernoulli Society, 2017, pp. 1326–1350  DOI: https://doi.org/10.1214/17-EJS1265
  * [9] Qianshun Cheng, Haiying Wang and Min Yang  “Information-based optimal subdata selection for big data logistic regression”  In _Journal of Statistical Planning and Inference_ 209 Elsevier, 2020, pp. 112–122  DOI: https://doi.org/10.1016/j.jspi.2020.03.004
  * [10] Stéphan Clémençon, Patrice Bertail and Emilie Chautru  “Scaling up M-estimation via sampling designs: The Horvitz-Thompson stochastic gradient descent”  In _2014 IEEE International Conference on Big Data (Big Data)_ , 2014, pp. 25–30  IEEE  DOI: https://doi.org/10.1109/BigData.2014.7004208
  * [11] S Cleveland and Ryan Hafen  “Divide and recombine (D& R): Data science for large complex data”  In _Statistical Analysis and Data Mining: The ASA Data Science Journal_ 7.6 Wiley Online Library, 2014, pp. 425–433  DOI: https://doi.org/10.1002/sam.11242
  * [12] Laura Deldossi and Chiara Tommasi  “Optimal design subsampling from Big Datasets”  In _Journal of Quality Technology_ 54.1 Taylor & Francis, 2022, pp. 93–101  DOI: https://doi.org/10.1080/00224065.2021.1889418
  * [13] Brian Donovan and Dan Work  “New York City Taxi Trip Data (2010-2013)”  University of Illinois at Urbana-Champaign, 2016  DOI: https://doi.org/10.13012/J8PN93H8
  * [14] Christopher C Drovandi, Christopher Holmes, James M McGree, Kerrie Mengersen, Sylvia Richardson and Elizabeth G Ryan  “Principles of Experimental Design for Big Data Analysis”  In _Statistical Science_ 32.3 Institute of Mathematical Statistics, 2017, pp. 385–404  DOI: https://doi.org/10.1214/16-STS604
  * [15] Ludwig Fahrmeir and Heinz Kaufmann  “Consistency and Asymptotic Normality of the Maximum Likelihood Estimator in Generalized Linear Models”  In _The Annals of Statistics_ 13.1 Institute of Mathematical Statistics, 1985, pp. 342–368  DOI: https://doi.org/10.1214/aos/1176346597
  * [16] Saptarshi Guha, Ryan Hafen, Jeremiah Rounds, Jin Xia, Jianfu Li, Bowei Xi and S Cleveland  “Large complex data: divide and recombine (D & R) with RHIPE”  In _Stat_ 1.1 Wiley Online Library, 2012, pp. 53–67  DOI: https://doi.org/10.1002/sta4.7
  * [17] Trevor Hastie and Robert Tibshirani  “Generalized Additive Models”  In _Statistical Science_ 1.3 Institute of Mathematical Statistics, 1986, pp. 297–310  URL: https://www.jstor.org/stable/2245459
  * [18] Bikram Karmakar and Indranil Mukhopadhyay  “Statistical Validity and Consistency of Big Data Analytics: A General Framework”  In _Statistics and Applications_ 18.2, 2020, pp. 369–381  URL: https://ssca.org.in/media/25_Vol._18_No._2_2020_SA_Indranil_Mukhopadhyay.pdf
  * [19] Ariel Kleiner, Ameet Talwalkar, Purnamrita Sarkar and I Jordan  “A scalable bootstrap for massive data”  In _Journal of the Royal Statistical Society: Series B: Statistical Methodology_ 76.4 JSTOR, 2014, pp. 795–816  DOI: https://doi.org/10.1111/rssb.12050
  * [20] Arvind Krishna, V Roshan Joseph, Shan Ba, William A Brenneman and William R Myers  “Robust experimental designs for model calibration”  In _Journal of Quality Technology_ 0.0 Taylor & Francis, 2021, pp. 1–12  DOI: https://doi.org/10.1080/00224065.2021.1930618
  * [21] JooChul Lee, Elizabeth D Schifano and HaiYing Wang  “Fast Optimal Subsampling Probability Approximation for Generalized Linear Models”  In _Econometrics and Statistics_ Elsevier, 2021  DOI: https://doi.org/10.1016/j.ecosta.2021.02.007
  * [22] Chengrui Li, Ying Hung and Minge Xie  “A sequential split-and-conquer approach for the analysis of big dependent data in computer experiments”  In _Canadian Journal of Statistics_ 48.4 Wiley Online Library, 2020, pp. 712–730  DOI: https://doi.org/10.1002/cjs.11559
  * [23] Nan Lin and Ruibin Xi  “Aggregated estimating equation estimation”  In _Statistics and Its Interface_ 4.1 International Press of Boston, 2011, pp. 73–83  DOI: https://doi.org/10.4310/SII.2011.v4.n1.a8
  * [24] Ping Ma, W Mahoney and Bin Yu  “A Statistical Perspective on Algorithmic Leveraging”  In _Journal of Machine Learning Research_ 16.27 JMLR. org, 2015, pp. 861–911  URL: http://jmlr.org/papers/v16/ma15a.html
  * [25] Ping Ma and Xiaoxiao Sun  “Leveraging for big data regression”  In _Wiley Interdisciplinary Reviews: Computational Statistics_ 7.1 Wiley Online Library, 2015, pp. 70–76  DOI: https://doi.org/10.1002/wics.1324
  * [26] Diana Martinez-Mosquera, Rosa Navarrete and Sergio Lujan-Mora  “Modeling and Management Big Data in Databases—A Systematic Literature Review”  In _Sustainability_ 12.2 Multidisciplinary Digital Publishing Institute, 2020, pp. 634  DOI: https://doi.org/10.3390/su12020634
  * [27] Cheng Meng, Rui Xie, Abhyuday Mandal, Xinlian Zhang, Wenxuan Zhong and Ping Ma  “LowCon: A Design-based Subsampling Approach in a Misspecified Linear Model”  In _Journal of Computational and Graphical Statistics_ 30.3 Taylor & Francis, 2021, pp. 694–708  DOI: https://doi.org/10.1080/10618600.2020.1844215
  * [28] John Ashworth Nelder and Robert WM Wedderburn  “Generalized Linear Models”  In _Journal of the Royal Statistical Society: Series A (General)_ 135.3 Wiley Online Library, 1972, pp. 370–384  DOI: https://doi.org/10.2307/2344614
  * [29] Trung Nguyen, Raymond G Gosine and Peter Warrian  “A Systematic Review of Big Data Analytics for Oil and Gas Industry 4.0”  In _IEEE access_ 8 IEEE, 2020, pp. 61183–61201  DOI: https://doi.org/10.1109/ACCESS.2020.2979678
  * [30] Bhatt Rajen and Dhall Abhinav  “Skin Segmentation”, UCI Machine Learning Repository, 2012  URL: https://archive.ics.uci.edu/ml/datasets/skin+segmentation
  * [31] Arshia Rehman, Saeeda Naz and Imran Razzak  “Leveraging big data analytics in healthcare enhancement: trends, challenges and opportunities”  In _Multimedia Systems_ Springer, 2021, pp. 1–33  DOI: https://doi.org/10.1007/s00530-020-00736-8
  * [32] D Schifano, Jing Wu, Chun Wang, Jun Yan and Ming-Hui Chen  “Online Updating of Statistical Inference in the Big Data Setting”  In _Technometrics_ 58.3 Taylor & Francis, 2016, pp. 393–403  DOI: https://doi.org/10.1080/00401706.2016.1142900
  * [33] Gideon Schwarz  “Estimating the Dimension of a Model”  In _The Annals of Statistics_ 6 Institute of Mathematical Statistics, 1978, pp. 461–464  URL: http://www.jstor.org/stable/2958889
  * [34] Chenlu Shi and Boxin Tang  “Model-Robust Subdata Selection for Big Data”  In _Journal of Statistical Theory and Practice_ 15.4 Springer, 2021, pp. 1–17  DOI: https://doi.org/10.1007/s42519-021-00217-9
  * [35] Aad W Van der Vaart  “Asymptotic Statistics” 3, Asymptotic Statistics  Cambridge university press, 2000
  * [36] Gregory Vaughan  “Efficient big data model selection with applications to fraud detection”  In _International Journal of Forecasting_ 36.3 Elsevier, 2020, pp. 1116–1127  DOI: https://doi.org/10.1016/j.ijforecast.2018.03.002
  * [37] Chun Wang, Ming-Hui Chen, Elizabeth Schifano, Jing Wu and Jun Yan  “Statistical methods and computing for big data”  In _Statistics and Its Interface_ 9.4 International Press of Boston, 2016, pp. 399–414  DOI: https://doi.org/10.4310/SII.2016.v9.n4.a1
  * [38] HaiYing Wang  “More Efficient Estimation for Logistic Regression with Optimal Subsamples”  In _Journal of Machine Learning Research_ 20.132, 2019, pp. 1–59  URL: http://jmlr.org/papers/v20/18-596.html
  * [39] HaiYing Wang, Min Yang and John Stufken  “Information-Based Optimal Subdata Selection for Big Data Linear Regression”  In _Journal of the American Statistical Association_ 114.525 Taylor & Francis, 2019, pp. 393–405  DOI: https://doi.org/10.1080/01621459.2017.1408468
  * [40] HaiYing Wang, Rong Zhu and Ping Ma  “Optimal Subsampling for Large Sample Logistic Regression”  In _Journal of the American Statistical Association_ 113.522 Taylor & Francis, 2018, pp. 829–844  DOI: https://doi.org/10.1080/01621459.2017.1292914
  * [41] Yishu Xue, HaiYing Wang, Jun Yan and Elizabeth D Schifano  “An online updating approach for testing the proportional hazards assumption with streams of survival data”  In _Biometrics_ 76.1 Wiley Online Library, 2020, pp. 171–182  DOI: https://doi.org/10.1111/biom.13137
  * [42] Ci Yang and Eric J Gonzales  “Modeling Taxi Trip Demand by Time of Day in New York City”  In _Transportation Research Record_ 2429.1 SAGE Publications Sage CA: Los Angeles, CA, 2014, pp. 110–120  DOI: https://doi.org/10.3141/2429-12
  * [43] Ci Yang and Eric J Gonzales  “Modeling Taxi Demand and Supply in New York City Using Large-Scale Taxi GPS Data”  In _Seeing Cities Through Big Data: Research, Methods and Applications in Urban Informatics_ Springer International Publishing, 2017, pp. 405–425  DOI: https://doi.org/10.1007/978-3-319-40902-3_22
  * [44] Yaqiong Yao and HaiYing Wang  “Optimal subsampling for softmax regression”  In _Statistical Papers_ 60.2 Springer, 2019, pp. 585–599  DOI: https://doi.org/10.1007/s00362-018-01068-6
  * [45] Yaqiong Yao and HaiYing Wang  “A Review on Optimal Subsampling Methods for Massive Datasets”  In _Journal of Data Science_ 19.1 School of Statistics, Renmin University of China, 2021, pp. 151–172  DOI: https://doi.org/10.6339/21-JDS999
  * [46] M Anil Yazici, Camille Kamga and Abhishek Singhal  “A big data driven model for taxi drivers’ airport pick-up decisions in new york city”  In _2013 IEEE International Conference on Big Data_ , 2013, pp. 37–44  IEEE  DOI: https://doi.org/10.1109/BigData.2013.6691775
  * [47] Jun Yu and HaiYing Wang  “Subdata selection algorithm for linear model discrimination”  In _Statistical Papers_ Springer, 2022, pp. 1–24  DOI: https://doi.org/10.1007/s00362-022-01299-8
  * [48] Yanxia Zhang and Yongheng Zhao  “Astronomy in the big data era”  In _Data Science Journal_ 14 Ubiquity Press, 2015  DOI: https://doi.org/10.5334/dsj-2015-011

##  Appendix A Proof of Theorem 5

###### Proof.

$\displaystyle\sum_{q=1}^{Q}\alpha_{q}\mbox{tr}{(\bm{V}_{q})}=\sum_{q=1}^{Q}\alpha_{q}\mbox{tr}\Big{(}\bm{J}^{-1}_{\bm{X}_{q}}{\bm{V}_{q}}_{c}\bm{J}^{-1}_{\bm{X}_{q}}\Big{)}$
$\displaystyle=$
$\displaystyle\frac{1}{N^{2}r}\sum_{q=1}^{Q}\alpha_{q}\sum_{i=1}^{N}\mbox{tr}\Big{[}\frac{1}{{\phi_{q}}_{i}}\\{y_{i}-\dot{\psi}(u(\hat{\bm{\theta}_{q}}^{T}_{MLE}\bm{x}_{qi}))\\}^{2}\bm{J}^{-1}_{\bm{X}_{q}}$
$\displaystyle\hskip
85.35826pt\dot{u}(\hat{\bm{\theta}_{q}}^{T}_{MLE}\bm{x}_{qi})\bm{x}_{qi}[\dot{u}(\hat{\bm{\theta}_{q}}_{MLE}^{T}\bm{x}_{qi})\bm{x}_{qi}]^{T}\bm{J}^{-1}_{\bm{X}_{q}}\Big{]}$
$\displaystyle=$
$\displaystyle\frac{1}{N^{2}r}\sum_{q=1}^{Q}\alpha_{q}\sum_{i=1}^{N}\Big{[}\frac{1}{{\phi_{q}}_{i}}\\{y_{i}-\dot{\psi}(u(\hat{\bm{\theta}_{q}}^{T}_{MLE}\bm{x}_{qi}))\\}^{2}||\bm{J}^{-1}_{\bm{X}_{q}}\dot{u}(\hat{\bm{\theta}_{q}}^{T}_{MLE}\bm{x}_{qi})\bm{x}_{qi}||^{2}\Big{]}$
$\displaystyle=$
$\displaystyle\frac{1}{N^{2}r}\sum_{q=1}^{Q}\alpha_{q}\Big{(}\sum_{i=1}^{N}{\phi_{q}}_{i}\Big{)}\sum_{i=1}^{N}\Big{[}\phi^{-1}_{qi}\\{y_{i}-\dot{\psi}(u(\hat{\bm{\theta}_{q}}^{T}_{MLE}\bm{x}_{qi}))\\}^{2}||\bm{J}^{-1}_{\bm{X}_{q}}\dot{u}(\hat{\bm{\theta}_{q}}^{T}_{MLE}\bm{x}_{qi})\bm{x}_{qi}||^{2}\Big{]}$
$\displaystyle=$
$\displaystyle\frac{\alpha_{1}}{N^{2}r}\Big{(}\sum_{i=1}^{N}{\phi_{1}}_{i}\Big{)}\sum_{i=1}^{N}\Big{[}\phi^{-1}_{1i}\\{y_{i}-\dot{\psi}(u(\hat{\bm{\theta}_{1}}^{T}_{MLE}\bm{x}_{1i}))\\}^{2}||\bm{J}^{-1}_{\bm{X}_{1}}\dot{u}(\hat{\bm{\theta}_{1}}^{T}_{MLE}\bm{x}_{1i})\bm{x}_{1i}||^{2}\Big{]}+\ldots+$
$\displaystyle\frac{\alpha_{Q}}{N^{2}r}\Big{(}\sum_{i=1}^{N}{\phi_{Q}}_{i}\Big{)}\sum_{i=1}^{N}\Big{[}\phi^{-1}_{Qi}\\{y_{i}-\dot{\psi}(u(\hat{\bm{\theta}_{Q}}^{T}_{MLE}\bm{x}_{Qi}))\\}^{2}||\bm{J}^{-1}_{\bm{X}_{Q}}\dot{u}(\hat{\bm{\theta}_{Q}}^{T}_{MLE}\bm{x}_{Qi})\bm{x}_{Qi}||^{2}\Big{]}$
$\displaystyle\geq$
$\displaystyle\frac{\alpha_{1}}{N^{2}r}\sum_{i=1}^{N}\\{y_{i}-\dot{\psi}(u(\hat{\bm{\theta}_{1}}^{T}_{MLE}\bm{x}_{1i}))\\}^{2}||\bm{J}^{-1}_{\bm{X}_{1}}\dot{u}(\hat{\bm{\theta}_{1}}^{T}_{MLE}\bm{x}_{1i})\bm{x}_{1i}||^{2}\Big{]}+\ldots+$
$\displaystyle\frac{\alpha_{Q}}{N^{2}r}\sum_{i=1}^{N}\\{y_{i}-\dot{\psi}(u(\hat{\bm{\theta}_{Q}}^{T}_{MLE}\bm{x}_{Qi}))\\}^{2}||\bm{J}^{-1}_{\bm{X}_{Q}}\dot{u}(\hat{\bm{\theta}_{Q}}^{T}_{MLE}\bm{x}_{Qi})\bm{x}_{Qi}||^{2}\Big{]}$
$\displaystyle=$
$\displaystyle\frac{1}{N^{2}r}\sum_{q=1}^{Q}\alpha_{q}\sum_{i=1}^{N}\\{y_{i}-\dot{\psi}(u(\hat{\bm{\theta}_{q}}^{T}_{MLE}\bm{x}_{qi}))\\}^{2}||\bm{J}^{-1}_{\bm{X}_{q}}\dot{u}(\hat{\bm{\theta}_{q}}^{T}_{MLE}\bm{x}_{qi})\bm{x}_{qi}||^{2}\Big{]}$
(13)

where the last inequality follows from the Cauchy-Schwarz inequality, and the
equality in it holds if and only if

${\phi_{q}}_{i}\propto|y_{i}-\dot{\psi}(u(\hat{\bm{\theta}_{q}}^{T}_{MLE}\bm{x}_{qi}))|\,||\bm{J}^{-1}_{\bm{X}_{q}}\dot{u}(\hat{\bm{\theta}_{q}}^{T}_{MLE}\bm{x}_{qi})\bm{x}_{qi}||.$

Here we define $0/0=0$, and this equivalent to removing data points with
$|y_{i}-\dot{\psi}(u(\hat{\bm{\theta}_{q}}^{T}_{MLE}\bm{x}_{qi}))|=0$ in the
expression of ${\bm{V}_{q}}_{c}$.

For the $q$-th model sub-sampling probabilities would be

${\phi_{q}}^{mMSE}_{i}=\frac{|y_{i}-\dot{\psi}(u(\hat{\bm{\theta}_{q}}^{T}_{MLE}\bm{x}_{qi}))|\,||\bm{J}^{-1}_{\bm{X}_{q}}\dot{u}(\hat{\bm{\theta}_{q}}^{T}_{MLE}\bm{x}_{qi})\bm{x}_{qi}||}{\sum_{j=1}^{N}|y_{j}-\dot{\psi}(u(\hat{\bm{\theta}_{q}}^{T}_{MLE}\bm{x}_{qj}))|\,||\bm{J}^{-1}_{\bm{X}_{q}}\dot{u}(\hat{\bm{\theta}_{q}}^{T}_{MLE}\bm{x}_{qj})\bm{x}_{qj}||},$
(14)

$i=1,\ldots,N$ and $\phi_{q}=\sum_{i=1}^{N}{\phi_{q}}_{i}^{mMSE}=1$.

If so, using the a priori probabilities $\alpha_{q}$ the model average optimal
sub-sampling probabilities is chosen such that

${\phi_{i}}^{mMSE}=\sum_{q=1}^{Q}\alpha_{q}{\phi_{q}}^{mMSE}_{i},$ (15)

$i=1,\ldots,N$, $\sum_{q=1}^{Q}\alpha_{q}=1$ then
$\sum_{q=1}^{Q}\alpha_{q}\mbox{tr}(\bm{V}_{q})$ attains its minimum.

Similarly optimal sub-sampling probabilities can be obtained for
$L$-optimality. ∎

##  Appendix B Algorithms

1

Stage 1

2 Random Sampling: Assign $\bm{\phi}=(\phi_{1},\ldots,\phi_{N})$ for $F_{N}$.
Here, $\bm{\phi}=\phi^{prop}$, where $\phi^{prop}_{i}=(2N_{0})^{-1}$ if the
$i$-th response element is zero, otherwise $\phi^{prop}_{i}=(2N_{1})^{-1}$,
and $N_{0}$ and $N_{1}$ are the number of elements in $\bm{y}$ when $y_{i}\in
0$ and $y_{i}\in 1$, respectively.

According to $\bm{\phi}$ draw a random sub-sample of size $r_{0}$,
$S_{r_{0}}=\\{\bm{x}^{r_{0}}_{l},y^{r_{0}}_{l},\phi^{r_{0}}_{l}\\}_{l=1}^{r_{0}}$.

Estimation: Based on $S_{r_{0}}$, using the Newton Raphson method find
$\tilde{\bm{\theta}}^{r_{0}}$ until $\tilde{\bm{\theta}}^{t+1}$ and
$\tilde{\bm{\theta}}^{t}$ are ‘close’ enough
($\tilde{\bm{\theta}}^{t+1}-\tilde{\bm{\theta}}^{t}<10^{-4}$),
$\displaystyle\tilde{\bm{\theta}}^{t+1}=\tilde{\bm{\theta}}^{t}+\Bigg{(}\sum_{l=1}^{r_{0}}\frac{\pi_{l}(\tilde{\bm{\theta}}^{t})(1-\pi_{l}(\tilde{\bm{\theta}}^{t}))\bm{x}^{r_{0}}_{l}(\bm{x}^{r_{0}}_{l})^{T}}{\phi^{r_{0}}_{l}}\Bigg{)}^{-1}\Bigg{(}\sum_{l=1}^{r_{0}}\frac{(y^{r_{0}}_{l}-\pi_{l}(\tilde{\bm{\theta}}^{t}))\bm{x}^{r_{0}}_{l}}{\phi^{r_{0}}_{l}}\Bigg{)},$
where $\mbox{logit}~{}\pi_{l}(\bm{\theta})=\bm{\theta}^{T}\bm{x}_{l}$.

3 Stage 2

4 Optimal sub-sampling probability: Using $\tilde{\bm{\theta}}^{r_{0}}$
estimate optimal sub-sampling probabilities,
$\displaystyle\phi^{mMSE}_{i}=\frac{|y_{i}-\pi_{i}(\tilde{\bm{\theta}}^{r_{0}})|\,||\tilde{\bm{J}}^{-1}_{\bm{X}}\bm{x}_{i}||}{\sum_{j=1}^{N}|y_{j}-\pi_{j}(\tilde{\bm{\theta}}^{r_{0}})|\,||\tilde{\bm{J}}^{-1}_{\bm{X}}\bm{x}_{j}||}\quad\mbox{or}\quad\phi^{mV_{c}}_{i}=\frac{|y_{i}-\pi_{i}(\tilde{\bm{\theta}}^{r_{0}})|\,||\bm{x}_{i}||}{\sum_{j=1}^{N}|y_{j}-\pi_{j}(\tilde{\bm{\theta}}^{r_{0}})|\,||\bm{x}_{j}||},$
where
$\tilde{\bm{J_{X}}}={(Nr_{0})}^{-1}\sum_{l=1}^{r_{0}}(\phi_{l}^{r_{0}})^{-1}\pi_{l}(\tilde{\bm{\theta}}^{r_{0}})(1-\pi_{l}(\tilde{\bm{\theta}}^{r_{0}})){\bm{x}_{l}}^{r_{0}}({\bm{x}_{l}}^{r_{0}})^{T}$
and $i=1,\ldots,N$.

5 Optimal sub-sampling and estimation: Using estimated $\bm{\phi}^{mMSE}$ or
$\bm{\phi}^{mV_{c}}$ draw a sub-sample $S_{r}$ completely at random (with
replacement) of size $r$ from $F_{N}$, such that
$S_{r}=\\{\bm{x}^{r}_{l},y^{r}_{l},\phi^{r}_{l}\\}_{l=1}^{r}$.

Form $S_{r_{0}+r}$ by combining $S_{r_{0}}$, $S_{r}$ and obtain
$\tilde{\bm{\theta}}$ until $\tilde{\bm{\theta}}^{t+1}$ and
$\tilde{\bm{\theta}}^{t}$ are ’close’ enough
($\tilde{\bm{\theta}}^{t+1}-\tilde{\bm{\theta}}^{t}<10^{-4}$):
$\displaystyle\tilde{\bm{\theta}}^{t+1}=\tilde{\bm{\theta}}^{t}+\Bigg{(}\sum_{k\in\\{r_{0},r\\}}\sum_{l=1}^{k}\frac{\pi_{l}(\tilde{\bm{\theta}}^{t})(1-\pi_{l}(\tilde{\bm{\theta}}^{t}))\bm{x}^{k}_{l}(\bm{x}^{k}_{l})^{T}}{\phi^{k}_{l}}\Bigg{)}^{-1}\Bigg{(}\sum_{k\in\\{r_{0},r\\}}\sum_{l=1}^{k}\frac{(y^{k}_{l}-\pi_{l}(\tilde{\bm{\theta}}^{t}))\bm{x}^{k}_{l}}{\phi^{k}_{l}}\Bigg{)}.$

Estimate the variance-covariance matrix of $\tilde{\bm{\theta}}$,
$\tilde{\bm{V}}$, by
$\tilde{\bm{J}}^{-1}_{\bm{X}}\tilde{\bm{V}}_{c}\tilde{\bm{J}}^{-1}_{\bm{X}}$,
where,  $\displaystyle\tilde{\bm{J_{X}}}=$
$\displaystyle\frac{1}{N(r_{0}+r)}\sum_{k\in\\{r_{0},r\\}}\sum_{l=1}^{k}\frac{\pi_{l}(\tilde{\bm{\theta}})(1-\pi_{l}(\tilde{\bm{\theta}})){\bm{x}^{k}_{l}}({\bm{x}^{k}_{l}})^{T}}{\phi^{k}_{l}}$
and
$\displaystyle\quad\tilde{\bm{V_{c}}}=\frac{1}{N^{2}(r_{0}+r)^{2}}\sum_{k\in\\{r_{0},r\\}}\sum_{l=1}^{k}\frac{(y^{k}_{l}-{\pi_{l}}(\tilde{\bm{\theta}}))^{2}{\bm{x}^{k}_{l}}({\bm{x}^{k}_{l}})^{T}}{({\phi^{k}_{l}})^{2}}.$

Output : $\tilde{\bm{\theta}},S_{r_{0}+r}$ and $\tilde{\bm{V}}$.

Algorithm 4 Two stage optimal sub-sampling algorithm for logistic regression.

1

Stage 1

2 Random Sampling: Assign $\bm{\phi}=(\phi_{1},\ldots,\phi_{N})$ for $F_{N}$.
Here, $\phi_{i}={N^{-1}}$.

According to $\bm{\phi}$ draw a random sub-sample of size $r_{0}$,
$S_{r_{0}}=\\{\bm{x}^{r_{0}}_{l},y^{r_{0}}_{l},\phi^{r_{0}}_{l}\\}_{l=1}^{r_{0}}=(\bm{X}^{r_{0}},\bm{y}^{r_{0}},\bm{\phi}^{r_{0}})$.

Estimation: Based on $S_{r_{0}}$, find:
$\displaystyle\tilde{\bm{\theta}}^{r_{0}}=\operatorname*{arg\,max}_{\bm{\theta}}~{}\log{L(\bm{\theta}|\bm{X}^{r_{0}},\bm{y}^{r_{0}},\bm{\phi}^{r_{0}})}\equiv\operatorname*{arg\,max}_{\bm{\theta}}~{}\frac{1}{r_{0}}\sum_{l=1}^{r_{0}}\Big{[}\frac{y^{r_{0}}_{l}\bm{\theta}^{T}\bm{x}^{r_{0}}_{l}-\lambda_{l}(\bm{\theta})-\log{(y^{r_{0}}_{l}!)}}{\phi^{r_{0}}_{l}}\Big{]},$
where $\log~{}{\lambda_{l}(\bm{\theta})}=\bm{\theta}^{T}\bm{x}_{l}.$

3 Stage 2

4 Optimal sub-sampling probability: Using $\tilde{\bm{\theta}}^{r_{0}}$
estimate optimal sub-sampling probabilities,
$\displaystyle\phi^{mMSE}_{i}=\frac{|y_{i}-\lambda_{i}(\tilde{\bm{\theta}}^{r_{0}})|\,||\tilde{\bm{J}}^{-1}_{\bm{X}}\bm{x}_{i}||}{\sum_{j=1}^{N}|y_{j}-\lambda_{j}(\tilde{\bm{\theta}}^{r_{0}})|\,||\tilde{\bm{J}}^{-1}_{\bm{X}}\bm{x}_{j}||}\quad\mbox{or}\quad\phi^{mV_{c}}_{i}=\frac{|y_{i}-\lambda_{i}(\tilde{\bm{\theta}}^{r_{0}})|\,||\bm{x}_{i}||}{\sum_{j=1}^{N}|y_{j}-\lambda_{j}(\tilde{\bm{\theta}}^{r_{0}})|\,||\bm{x}_{j}||},$
where
$\tilde{\bm{J_{X}}}={(Nr_{0})}^{-1}\sum_{l=1}^{r_{0}}(\phi_{l}^{r_{0}})^{-1}\lambda_{l}(\tilde{\bm{\theta}}^{r_{0}}){\bm{x}_{l}}^{r_{0}}({\bm{x}_{l}}^{r_{0}})^{T}$
and $i=1,\ldots,N$.

5 Optimal sub-sampling and estimation: Using estimated $\bm{\phi}^{mMSE}$ or
$\bm{\phi}^{mV_{c}}$ draw a sub-sample $S_{r}$ completely at random (with
replacement) of size $r$ from $F_{N}$, such that
$S_{r}=\\{\bm{x}^{r}_{l},y^{r}_{l},\phi^{r}_{l}\\}_{l=1}^{r}=(\bm{X}^{r},\bm{y}^{r},\bm{\phi}^{r})$.

Form $S_{r_{0}+r}$ by combining $S_{r_{0}}$, $S_{r}$ and obtain:
$\displaystyle\tilde{\bm{\theta}}=$
$\displaystyle\operatorname*{arg\,max}_{\bm{\theta}}~{}\log{L(\bm{\theta}|\bm{X}^{r_{0}},\bm{y}^{r_{0}},\bm{\phi}^{r_{0}},\bm{X}^{r},\bm{y}^{r},\bm{\phi}^{r})}$
$\displaystyle\equiv$
$\displaystyle\operatorname*{arg\,max}_{\bm{\theta}}~{}\frac{1}{r_{0}+r}\sum_{k\in\\{r_{0},r\\}}\sum_{l=1}^{k}\Big{[}\frac{y^{k}_{l}\bm{\theta}^{T}\bm{x}^{k}_{l}-\lambda_{l}(\bm{\theta})-\log{(y^{k}_{l}!)}}{\phi^{k}_{l}}\Big{]}.$

Estimate the variance-covariance matrix of $\tilde{\bm{\theta}}$,
$\tilde{\bm{V}}$, by
$\tilde{\bm{J}}^{-1}_{\bm{X}}\tilde{\bm{V}}_{c}\tilde{\bm{J}}^{-1}_{\bm{X}}$,
where,  $\displaystyle\tilde{\bm{J_{X}}}=$
$\displaystyle\frac{1}{N(r_{0}+r)}\sum_{k\in\\{r_{0},r\\}}\sum_{l=1}^{k}\frac{\lambda_{l}(\tilde{\bm{\theta}}){\bm{x}^{k}_{l}}({\bm{x}^{k}_{l}})^{T}}{\phi^{k}_{l}}$
and
$\displaystyle\quad\tilde{\bm{V_{c}}}=\frac{1}{N^{2}(r_{0}+r)^{2}}\sum_{k\in\\{r_{0},r\\}}\sum_{l=1}^{k}\frac{(y^{k}_{l}-{\lambda_{l}}(\tilde{\bm{\theta}}))^{2}{\bm{x}^{k}_{l}}({\bm{x}^{k}_{l}})^{T}}{({\phi^{k}_{l}})^{2}}.$

Output : $\tilde{\bm{\theta}},S_{r_{0}+r}$ and $\tilde{\bm{V}}$.

Algorithm 5 Two stage optimal sub-sampling algorithm for Poisson regression.

1

Stage 1

2 Random Sampling: Assign $\bm{\phi}=(\phi_{1},\ldots,\phi_{N})$ for $F_{N}$.
Here, $\bm{\phi}=\phi^{prop}$, where $\phi^{prop}_{i}=(2N_{0})^{-1}$ if the
$i$-th response element is $0$, otherwise $\phi^{prop}_{i}=(2N_{1})^{-1}$, and
$N_{0}$ and $N_{1}$ are the number of elements in $\bm{y}$ when $y_{i}\in 0$
and $y_{i}\in 1$, respectively.

According to $\bm{\phi}$ draw random sub-samples of size $r_{0}$, such that
$S_{r_{0}}=\\{h_{q}(\bm{{x_{0}}}^{r_{0}}_{l}),y^{r_{0}}_{l},\phi^{r_{0}}_{l}\\}_{l=1}^{r_{0}}$
for $q=1,\ldots,Q$.

Estimation: For $q=1,\ldots,Q$ and $S_{r_{0}}$, using the Newton Raphson
method find $\tilde{\bm{\theta}_{q}}^{r_{0}}$ until
$\tilde{\bm{\theta}_{q}}^{t+1}$ and $\tilde{\bm{\theta}_{q}}^{t}$ are ‘close’
enough ($\tilde{\bm{\theta}_{q}}^{t+1}-\tilde{\bm{\theta}_{q}}^{t}<10^{-4}$),
$\displaystyle\tilde{\bm{\theta}_{q}}^{t+1}=\tilde{\bm{\theta}_{q}}^{t}+\Bigg{(}\sum_{l=1}^{r_{0}}\frac{\pi_{l}(\tilde{\bm{\theta}_{q}}^{t})(1-\pi_{l}(\tilde{\bm{\theta}_{q}}^{t}))\bm{x}_{ql}(\bm{x}_{ql})^{T}}{\phi_{l}}\Bigg{)}^{-1}\Bigg{(}\sum_{l=1}^{r_{0}}\frac{(y_{l}-\pi_{l}(\tilde{\bm{\theta}_{q}}^{t}))\bm{x}_{ql}}{\phi_{l}}\Bigg{)},$
where
$\mbox{logit}~{}\pi_{l}(\bm{\theta}_{q})=\bm{\theta}^{T}_{q}\bm{x}_{ql}$.

3 Stage 2

4 Optimal sub-sampling probability: Based on $\tilde{\bm{\theta}_{q}}^{r_{0}}$
and $\alpha_{q}$ (for $q=1,\ldots,Q$) estimate optimal sub-sampling
probabilities,
$\phi^{mMSE}_{i}=\sum_{q=1}^{Q}\frac{\alpha_{q}|y_{i}-\pi_{i}(\tilde{\bm{\theta}_{q}}^{r_{0}})|\,||\tilde{\bm{J}}^{-1}_{\bm{X}_{q}}\bm{x}_{qi}||}{\sum_{j=1}^{N}|y_{j}-\pi_{j}(\tilde{\bm{\theta}_{q}}^{r_{0}})|\,||\tilde{\bm{J}}^{-1}_{\bm{X}_{q}}\bm{x}_{qj}||}\quad\mbox{or}\quad\phi^{mV_{c}}_{i}=\sum_{q=1}^{Q}\frac{\alpha_{q}|y_{i}-\pi_{i}(\tilde{\bm{\theta}_{q}}^{r_{0}})|\,||\bm{x}_{qi}||}{\sum_{j=1}^{N}|y_{j}-\pi_{j}(\tilde{\bm{\theta}_{q}}^{r_{0}})|\,||\bm{x}_{qj}||},$
where
$\tilde{\bm{J}}_{\bm{X}_{q}}={(Nr_{0})}^{-1}\sum_{l=1}^{r_{0}}(\phi_{l}^{r_{0}})^{-1}\pi_{l}(\tilde{\bm{\theta}_{q}}^{r_{0}})(1-\pi_{l}(\tilde{\bm{\theta}_{q}}^{r_{0}})){\bm{x}^{r_{0}}_{ql}}(\bm{x}^{r_{0}}_{ql})^{T}$
and $i=1,\ldots,N$.

5 Optimal sub-sampling and estimation: Based on $\bm{\phi}^{mMSE}$ or
$\bm{\phi}^{mV_{c}}$ draw a sub-sample of size $r$ from $F_{N}$, such that
$S_{r}=\\{h_{q}(\bm{{x_{0}}}^{r}_{l}),y^{r}_{l},\phi^{r}_{l}\\}_{l=1}^{r}$ for
$q=1,\ldots,Q$.

Combine $S_{r_{0}},S_{r}$ and form $S_{(r_{0}+r)}$, and obtain
$\tilde{\bm{\theta}_{q}}$ until $\tilde{\bm{\theta}_{q}}^{t+1}$ and
$\tilde{\bm{\theta}_{q}}^{t}$ are ’close’ enough
($\tilde{\bm{\theta}_{q}}^{t+1}-\tilde{\bm{\theta}_{q}}^{t}<10^{-4}$):
$\displaystyle\tilde{\bm{\theta}_{q}}^{t+1}=\tilde{\bm{\theta}_{q}}^{t}+\Bigg{[}\sum_{k\in\\{r_{0},r\\}}\sum_{l=1}^{k}\frac{\pi_{l}(\tilde{\bm{\theta}_{q}}^{t})(1-\pi_{l}(\tilde{\bm{\theta}_{q}}^{t}))\bm{x}^{k}_{ql}(\bm{x}^{k}_{ql})^{T}}{\phi^{k}_{l}}\Bigg{]}^{-1}\Bigg{[}\sum_{k\in\\{r_{0},r\\}}\sum_{l=1}^{k}\frac{(y^{k}_{l}-\pi_{l}(\tilde{\bm{\theta}_{q}}^{t}))\bm{x}^{k}_{ql}}{\phi^{k}_{l}}\Bigg{]}.$

Estimate the variance-covariance matrix of $\tilde{\bm{\theta}}_{q}$,
$\tilde{\bm{V}}_{q}$, by
$\tilde{\bm{J}}^{-1}_{\bm{X}_{q}}\tilde{\bm{V}}_{q_{c}}\tilde{\bm{J}}^{-1}_{\bm{X}_{q}}$,
where,  $\displaystyle\tilde{\bm{J}}_{\bm{X}_{q}}=$
$\displaystyle\frac{1}{N(r_{0}+r)}\sum_{k\in\\{r_{0},r\\}}\sum_{l=1}^{k}\frac{\pi_{l}(\tilde{\bm{\theta}_{q}})(1-\pi_{l}(\tilde{\bm{\theta}_{q}}))\bm{x}^{k}_{ql}(\bm{x}^{k}_{ql})^{T}}{\phi^{k}_{l}}$
and
$\displaystyle\quad\tilde{\bm{V}}_{q_{c}}=\frac{1}{N^{2}(r_{0}+r)^{2}}\sum_{k\in\\{r_{0},r\\}}\sum_{l=1}^{k}\frac{(y^{k}_{l}-{\pi_{l}}(\tilde{\bm{\theta}_{q}}))^{2}\bm{x}^{k}_{ql}(\bm{x}^{k}_{ql})^{T}}{({\phi^{k}_{l}})^{2}}.$

Output : $\tilde{\bm{\theta}}_{q},S_{(r_{0}+r)}$ and $\tilde{\bm{V}}_{q}$ for
$q=1,\ldots,Q$.

Algorithm 6 Two stage model robust optimal sub-sampling algorithm for logistic
regression.

1

Stage 1

2 Random Sampling : Assign $\bm{\phi}=(\phi_{1},\ldots,\phi_{N})$ for $F_{N}$.
Here, $\phi_{i}={N^{-1}}$.

According to $\bm{\phi}$ draw random sub-samples of size $r_{0}$, such that
$S_{r_{0}}=\\{h_{q}({\bm{x}_{0}}^{r_{0}}_{l}),y^{r_{0}}_{l},\phi^{r_{0}}_{l}\\}_{l=1}^{r_{0}}=(h_{q}(\bm{X}_{0}^{r_{0}}),\bm{y}^{r_{0}},\bm{\phi}^{r_{0}})$
for $q=1,\ldots,Q$.

Estimation: For $q=1,\ldots,Q$ and $S_{r_{0}}$ find:
$\displaystyle\tilde{\bm{\theta}_{q}}^{r_{0}}=\operatorname*{arg\,max}_{\bm{\theta}_{q}}~{}\log{L(\bm{\theta}_{q}|\bm{X}^{r_{0}}_{q},\bm{y}^{r_{0}},\bm{\phi}^{r_{0}})}\equiv\operatorname*{arg\,max}_{\bm{\theta}_{q}}~{}\frac{1}{r_{0}}\sum_{l=1}^{r_{0}}\Big{[}\frac{y^{r_{0}}_{l}\bm{\theta}^{T}_{q}\bm{x}^{r_{0}}_{ql}-\lambda_{l}(\bm{\theta}_{q})-\log{(y^{r_{0}}_{l}!)}}{\phi^{r_{0}}_{l}}\Big{]},$
where $\log{\lambda_{l}(\bm{\theta}_{q})}=\bm{\theta}^{T}_{q}\bm{x}_{ql}.$

3 Stage 2

4 Optimal sub-sampling probability: For $q=1,\ldots,Q$, using
$\tilde{\bm{\theta}_{q}^{r_{0}}}$ estimate optimal sub-sampling probabilities,
$\displaystyle\phi^{mMSE}_{i}=\sum_{q=1}^{Q}\frac{\alpha_{q}|y_{i}-\lambda_{i}(\tilde{\bm{\theta}_{q}}^{r_{0}})|\,||\tilde{\bm{J}}^{-1}_{\bm{X}_{q}}\bm{x}_{qi}||}{\sum_{j=1}^{N}|y_{j}-\lambda_{j}(\tilde{\bm{\theta}_{q}}^{r_{0}})|\,||\tilde{\bm{J}}^{-1}_{\bm{X}_{q}}\bm{x}_{qj}||}\quad\mbox{or}\quad\phi^{mV_{c}}_{i}=\sum_{q=1}^{Q}\frac{\alpha_{q}|y_{i}-\lambda_{i}(\tilde{\bm{\theta}_{q}}^{r_{0}})|\,||\bm{x}_{qi}||}{\sum_{j=1}^{N}|y_{j}-\lambda_{j}(\tilde{\bm{\theta}_{q}}^{r_{0}})|\,||\bm{x}_{qj}||},$
where
$\tilde{\bm{J}}_{\bm{X}_{q}}={(Nr_{0})}^{-1}\sum_{l=1}^{r_{0}}(\phi_{l}^{r_{0}})^{-1}\lambda_{l}(\tilde{\bm{\theta}_{q}}^{r_{0}})\bm{x}^{r_{0}}_{ql}(\bm{x}^{r_{0}}_{ql})^{T}$
and $i=1,\ldots,N$.

5 Optimal sub-sampling and estimation: Based on $\bm{\phi}^{mMSE}$ or
$\bm{\phi}^{mV_{c}}$ draw a sub-sample of size $r$ from $F_{N}$, such that
$S_{r}=\\{h_{q}({\bm{x}_{0}}^{r}_{l}),y^{r}_{l},\phi^{r}_{l}\\}_{l=1}^{r}=(h_{q}(\bm{X}_{0}^{r}),\bm{y}^{r}_{q},\bm{\phi}^{r}_{q})$
for $q=1,\ldots,Q$.

Combine $S_{r_{0}}$, $S_{r}$ and form $S_{(r_{0}+r)}$, and obtain:
$\displaystyle\tilde{\bm{\theta}}_{q}=$
$\displaystyle\operatorname*{arg\,max}_{\bm{\theta}_{q}}~{}\log{L(\bm{\theta}_{q}|\bm{X}^{r_{0}}_{q},\bm{y}^{r_{0}}_{q},\bm{\phi}^{r_{0}}_{q},\bm{X}^{r}_{q},\bm{y}^{r}_{q},\bm{\phi}^{r}_{q})}$
$\displaystyle\equiv$
$\displaystyle\operatorname*{arg\,max}_{\bm{\theta}_{q}}~{}\frac{1}{r_{0}+r}\sum_{k\in\\{r_{0},r\\}}\sum_{l=1}^{k}\Big{[}\frac{y^{k}_{l}\bm{\theta}^{T}_{q}\bm{x}^{k}_{ql}-\lambda_{l}(\bm{\theta}_{q})-\log{(y^{k}_{l}!)}}{\phi^{k}_{l}}\Big{]}.$

Estimate the variance-covariance matrix of $\tilde{\bm{\theta}}_{q}$,
$\tilde{\bm{V}}_{q}$, by
$\tilde{\bm{J}}^{-1}_{\bm{X}_{q}}\tilde{\bm{V}}_{q_{c}}\tilde{\bm{J}}^{-1}_{\bm{X}_{q}}$,
here,  $\displaystyle\tilde{\bm{J}}_{\bm{X}_{q}}=$
$\displaystyle\frac{1}{N(r_{0}+r)}\sum_{k\in\\{r_{0},r\\}}\sum_{l=1}^{k}\frac{\lambda_{l}(\tilde{\bm{\theta}_{q}})\bm{x}^{k}_{ql}(\bm{x}^{k}_{ql})^{T}}{\phi^{k}_{l}}$
and
$\displaystyle\quad\tilde{\bm{V}}_{q_{c}}=\frac{1}{N^{2}(r_{0}+r)^{2}}\sum_{k\in\\{r_{0},r\\}}\sum_{l=1}^{k}\frac{(y^{k}_{l}-{\lambda_{l}}(\tilde{\bm{\theta}_{q}}))^{2}\bm{x}^{k}_{ql}(\bm{x}^{k}_{ql})^{T}}{({\phi^{k}_{l}})^{2}}.$

Output : $\tilde{\bm{\theta}}_{q},S_{(r_{0}+r)}$ and $\tilde{\bm{V}}_{q}$ for
$q=1,\ldots,Q$.

Algorithm 7 Two stage model robust optimal sub-sampling algorithm for Poisson
regression.

##  Appendix C Code materials for the simulation setup and real world
applications

  1. 1.

Logistic Regression

    1. (a)

Two covariate model with covariate data generated by the Exponential
distribution.

    2. (b)

Two covariate model with covariate data generated by the Normal distribution.

    3. (c)

Skin segmentation data.

  2. 2.

Poisson Regression

    1. (a)

Two covariate model with covariate data generated by the Normal distribution.

    2. (b)

Two covariate model with covariate data generated by the Uniform distribution.

    3. (c)

New York City taxi fare data.

##  Appendix D Figures

Figure 5: (a) Model information (larger is better) and (b) SMSE (smaller is
better) for the sub-sampling methods for the logistic regression model under
$mMSE$. Covariate data were generated from the Exponential distribution.
Figure 6: (a) Model information (larger is better) and (b) SMSE (smaller is
better) for the sub-sampling methods for the Poisson regression model under
$mMSE$. Covariate data were generated from the Uniform distribution.