# Asymptotic Results for Penalized Quasi-Likelihood Estimation in Generalized
Linear Mixed Models

Xu Ning Corresponding author: Email<EMAIL_ADDRESS>Address: Research
School of Finance, Actuarial Studies and Statistics, The Australian National
University, Canberra, ACT, 2601, Australia Research School of Finance,
Actuarial Studies and Statistics, The Australian National University  Francis
K.C. Hui  Research School of Finance, Actuarial Studies and Statistics, The
Australian National University  A. H. Welsh  Research School of Finance,
Actuarial Studies and Statistics, The Australian National University

###### Abstract

Generalized Linear Mixed Models (GLMMs) are widely used for analysing
clustered data. One well-established method of overcoming the integral in the
marginal likelihood function for GLMMs is penalized quasi-likelihood (PQL)
estimation, although to date there are few asymptotic distribution results
relating to PQL estimation for GLMMs in the literature. In this paper, we
establish large sample results for PQL estimators of the parameters and random
effects in independent-cluster GLMMs, when both the number of clusters and the
cluster sizes go to infinity. This is done under two distinct regimes:
conditional on the random effects (essentially treating them as fixed effects)
and unconditionally (treating the random effects as random). Under the
conditional regime, we show the PQL estimators are asymptotically normal
around the true fixed and random effects. Unconditionally, we prove that while
the estimator of the fixed effects is asymptotically normally distributed, the
correct asymptotic distribution of the so-called prediction gap of the random
effects may in fact be a normal scale-mixture distribution under certain
relative rates of growth. A simulation study is used to verify the finite
sample performance of our theoretical results.

Keywords  Asymptotic independence, Clustered data, Large sample distribution,
Longitudinal data, Prediction

##  1 Introduction

Generalized linear mixed models (GLMMs) are widely used in statistics to model
relationships in clustered and correlated data [52]. As the marginal
likelihood function of GLMMs, except for normally distributed responses with
the identity link, contains an intractable integral, many methods have been
developed to estimate and perform inference for the parameters in a
computationally efficient manner. These include the Laplace approximation,
Gauss-Hermite quadrature, and variational approximations, among others [52,
56, 36]. A connected and well-established approach is penalized Quasi-
Likelihood (PQL) estimation [34]. As one of the first methods to circumvent
the intractable integral, PQL estimation has seen a resurgence in modern
statistics as a very computationally efficient method for high-dimensional
multivariate GLMMs [42, 49, e.g.,]. However, despite its long history, large
sample distributional results for PQL estimation in mixed models are scarce.

The most often studied asymptotic results for maximum likelihood estimators of
GLMMs are based on increasing the number of clusters while keeping the size of
each cluster fixed or bounded [52, 53]. Asymptotic results when both the
cluster size and number of clusters grow are less developed, although some
results for the maximum likelihood estimator as well as the empirical best
linear unbiased predictor (EBLUP) for the linear mixed model (LMM) have been
developed; see [50, 51] and references therein. Recently, [47] proved an
asymptotic normality result for a maximum quasi-likelihood estimator of the
fixed parameters, which is different from the PQL estimator, for independent-
cluster GLMMs.

This work is distinct from the above results: compared to [50, 51] we consider
a more general random effects structure that permits random slopes in GLMMs.
Meanwhile, [47] considered GLMMs but not the case when cluster sizes grow
faster than the number of clusters; nor did they present results for
predictors of random effects, both of which are considered in this article.
Furthermore, we establish results for the prediction gap in GLMMs, which are
new to the literature and allow unconditional inference to be performed for
the random effects [48, 58, noting unconditional inference for random effects
in LMMs has been considered previously in a very different way through the
unconditional mean squared error of prediction,]. Note for the PQL estimator
specifically, [60], [43] and [42] demonstrated estimation consistency under
increasing cluster size and number of clusters, but did not develop any large
sample distributional results.

It is important to remark that when cluster sizes do not increase, PQL is
known to be asymptotically biased [35, e.g.,]. As such, increasing both the
number of clusters and cluster size is a necessary condition for the PQL
estimator to be consistent. Indeed, increasing number of clusters and cluster
size is necessary for the consistency of other estimators such as the Laplace
estimator [54, 42, 55]. With this in mind, we develop our large sample
distributional results under this setting, with the precise rates of growth to
be formalised later. We note this asymptotic framework is relevant for many
applications with large cluster sizes e.g., educational studies with large
numbers of students (units) grouped within schools (clusters), and medical
studies with large groups (clusters) of patients (units) treated at different
hospitals.

We derive our asymptotic results for the PQL estimator under two distinct
sampling regimes. In the first, we condition on the random effects, i.e. treat
them as fixed effects, although we will still refer to them as random effects
for consistency. In the second, unconditional regime, we allow the random
effects to be random. Conditional inference is appropriate when hypothetical
resampling occurs within the same observed clusters, while unconditional
inference may be more appropriate when (new) clusters are sampled from some
population. Importantly, we demonstrate the asymptotic distributional results
for the two regimes differ markedly. Conditional on the random effects we show
the PQL estimator is asymptotically normally distributed around the true
parameter values, with a convergence rate of ${N}^{1/2}$ (square root of the
total number of observations) for the fixed effects and $n_{i}^{1/2}$ (square
root of the cluster size of the $i$th cluster) for the random effects (which
are now also fixed parameters). We find that when a variable is included as
both a fixed and random effect covariate, the PQL estimator is asymptotically
normally distributed around a sum-to-zero reparametrized version of the
estimand. Unconditionally, we demonstrate the asymptotic normality of the PQL
estimator for the fixed effects around the true values, but with a slower
convergence rate of ${m}^{1/2}$ (square root of the number of clusters).
Furthermore, we demonstrate that the “prediction gap” i.e., the difference
between the PQL estimator of and the true random effect, is not in general
asymptotically normally distributed; instead, it follows a normal scale-
mixture when $m$ grows faster than $n_{i}$.

Our results have important implications for inference in GLMMs. There is a
choice of whether conditional or unconditional inference is desired, with
different asymptotic distributions needing to be applied in each case. Also,
the potential asymptotic non-normality of the prediction gap has consequences
in practice, since normality is often assumed when constructing prediction
intervals for the random effects in GLMMs [32, 36]. The theoretical results in
this paper offer an important step towards more formal, rigorous asymptotic
inference using PQL estimation (and perhaps other similar estimators) for
GLMMs.

The structure of the article is as follows. In Section 2, we introduce GLMMs
and PQL estimation. In Sections 3 and 4, we present and develop our asymptotic
framework and results for the conditional and unconditional regimes. In
Section 5, we present results from a simulation study which empirically verify
our large sample developments. Finally, in Section 6 we discuss the
implications of our results.

##  2 Generalized Linear Mixed Models

We study the independent-cluster generalized linear mixed model defined as
follows. Let $y_{ij}$ denote the $j$th measurement of cluster $i$,
$\bm{x}_{ij}$ denote a vector of $p_{f}$ fixed effect covariates, and
$\bm{z}_{ij}$ denote a vector of $p_{r}$ random effect covariates, for
$j=1,\ldots,n_{i}$, and $i=1,\ldots,m$. Let $N=\sum_{i=1}^{m}n_{i}$,
$n=m^{-1}N$, $n_{L}=\min_{1\leq i\leq m}n_{i}$, and $n_{U}=\max_{1\leq i\leq
m}n_{i}$. The $m$ clusters are independent of each other. Conditional on a
$p_{r}$-vector of random effects $\dot{\bm{b}}_{i}$, where the dot above any
quantity is used to denote its true value (or that it is evaluated at the true
parameter values), the responses $y_{ij}$ from cluster $i$ are assumed to be
independent observations from the exponential family with mean
$\dot{\mu}_{ij}$ and dispersion parameter $\dot{\phi}$. That is,
$f(y_{ij}|\dot{\bm{\beta}},\dot{\bm{b}}_{i},\dot{\phi})=\exp\,[\dot{\phi}^{-1}\\{y_{ij}\dot{\vartheta}_{ij}-a(\dot{\vartheta}_{ij})\\}+c(y_{ij},\dot{\phi})]$,
for known functions $a(\cdot)$, $c(\cdot)$, and $g(\cdot)$ satisfying
$g(\dot{\mu}_{ij})=g\\{a^{\prime}(\dot{\vartheta}_{ij})\\}=\dot{\eta}_{ij}=\bm{x}^{\top}_{ij}\dot{\bm{\beta}}+\bm{z}^{\top}_{ij}\dot{\bm{b}}_{i}$,
where $\dot{\bm{\beta}}$ denotes a $p_{f}$-vector of true fixed effect
coefficients, and $\dot{\eta}_{ij}$ the corresponding true linear predictor.
For ease of development, we assume that the canonical link is used, so that
$\dot{\vartheta}=\dot{\eta}$. Commonly used distributions within the
exponential family include the normal, Poisson, binomial and gamma
distributions. The true realised random effects $\dot{\bm{b}}_{i}$ are
independently and identically distributed (i.i.d.), drawn from a multivariate
normal distribution with zero mean vector and unstructured $p_{r}\times p_{r}$
random effects covariance matrix $\dot{\bm{G}}$. That is,
$\dot{\bm{b}}_{i}\overset{i.i.d.}{\sim}N(\bm{0},\dot{\bm{G}})$.

Write $\bm{X}_{i}=[\bm{x}_{i1},\ldots,\bm{x}_{in_{i}}]^{\top}$, and
$\bm{Z}_{i}=[\bm{z}_{i1},\ldots,\bm{z}_{in_{i}}]^{\top}$, so we can
concatenate the means across the measurements for each cluster to obtain
$g(\dot{\bm{\mu}}_{i})=\bm{X}_{i}\dot{\bm{\beta}}+\bm{Z}_{i}\dot{\bm{b}}_{i}$
for $\dot{\bm{\mu}}_{i}=(\dot{\mu}_{i1},\ldots,\dot{\mu}_{in_{i}})^{\top}$,
where $g(\dot{\bm{\mu}}_{i})$ denotes applying the link function $g(\cdot)$ to
$\dot{\bm{\mu}}_{i}$ component-wise. We can further concatenate across
clusters and write
$g(\dot{\bm{\mu}})=\bm{X}\dot{\bm{\beta}}+\bm{Z}\dot{\bm{b}}$, with
$\dot{\bm{\mu}}=(\dot{\bm{\mu}}_{1}^{\top},\ldots,\dot{\bm{\mu}}_{m}^{\top})^{\top}$,
$\bm{X}=[\bm{X}^{\top}_{1},\ldots,\bm{X}^{\top}_{m}]^{\top}$,
$\bm{Z}=\text{bdiag}(\bm{Z}_{1},\ldots\bm{Z}_{m})$, and
$\dot{\bm{b}}=(\dot{\bm{b}}_{1}^{\top},\ldots,\dot{\bm{b}}_{m}^{\top})^{\top}$.
Here, $\text{bdiag}()$ is the block-diagonal matrix operator, $\bm{X}$ is of
dimension $N\times p_{f}$, and $\bm{Z}$ is an $N\times mp_{r}$ sparse block-
diagonal matrix, with at most $p_{r}$ non-zero components per row, and at most
$n_{U}$ non-zero components per column.

Let $\bm{y}_{i}=(y_{11},\ldots,y_{1n_{i}})^{\top}$ and
$\bm{y}=(\bm{y}_{1}^{\top},\ldots,\bm{y}_{m}^{\top})^{\top}$. Then the
marginal log-likelihood function for the independent-cluster GLMM is given by

$\displaystyle\ln f(\bm{y}|\bm{\beta},\phi,\bm{G})$
$\displaystyle=\sum_{i=1}^{m}\ln\int\left(\prod_{j=1}^{n_{i}}f(y_{ij}|\bm{\beta},\bm{b}_{i},\phi)\right)f(\bm{b}_{i}|\bm{G})\,d\bm{b}_{i}.$
(1)

The above integral is not available analytically except in the special case of
a normal response with an identity link function. Let
$\bm{\theta}=(\bm{\beta}^{\top},\bm{b}^{\top})^{\top}$ denote the full vector
of fixed and random effects. Then for a given $\bm{G}$ and $\phi$, the PQL
objective function for an independent-cluster GLMM is defined as

$\displaystyle Q(\bm{\theta})=\sum_{i=1}^{m}\sum_{j=1}^{n_{i}}\ln
f(y_{ij}|\bm{\beta},\bm{b}_{i},\phi)-\frac{1}{2}\sum_{i=1}^{m}\bm{b}_{i}^{\top}\bm{G}^{-1}\bm{b}_{i},$
(2)

and the PQL estimator is defined as
$\hat{\bm{\theta}}=\underset{\bm{\theta}}{\arg\max}\,Q(\bm{\theta})$. As there
are no integrals in (2), the computational cost of PQL estimation is low
relative to standard maximum likelihood estimation [34]. Note for normal
linear mixed models, the integral in the likelihood already possesses an
analytical solution when an identity link is used, and PQL estimation is
equivalent to the mixed model equations of [41] assuming the error variance is
known.

The PQL procedure provides explicit estimators of both the fixed and random
effects. The latter is practically useful since the random effects play an
important implicit role in fitting and using the GLMM. For instance, the
realised values of the random effects (or functions thereof) are often
important in prediction problems such as small-area estimation [45, 57], while
the empirical distribution of the random effects estimators is often examined
in model diagnostics [44]. On the other hand, (2) alone does not incorporate
estimation of the random effects covariance matrix. From a theoretical
standpoint, existing papers on large sample theory for PQL and related
objective functions have assumed $\dot{\bm{G}}$ is known for the purposes of
asymptotic development [60, 53, e.g.,]. Practically speaking, several
approaches have been suggested to estimate $\dot{\bm{G}}$ when applying PQL,
for example by using a working LMM [34], the Laplace objective function [43],
or simply the sample covariance matrix of the estimated random effects [46].
Indeed, [46] and [43] demonstrated that the sample covariance of the estimated
random effects is a consistent estimator of $\dot{\bm{G}}$ under suitable
regularity conditions.

In this article, we set $\bm{G}=\hat{\bm{G}}$ in (2), where $\hat{\bm{G}}$ is
a symmetric positive definite matrix that is either non-stochastic or its
inverse $\hat{\bm{G}}^{-1}$ is stochastically bounded. Importantly, our large
sample developments do not require $\hat{\bm{G}}$ to necessarily be a
consistent estimator of the true random effects covariance matrix
$\dot{\bm{G}}$. For example, while we can use the estimators of $\dot{\bm{G}}$
mentioned above, our theory also permits setting $\hat{\bm{G}}$ to some fixed
matrix e.g., the identity matrix, say. Intuitively, this is because we develop
our large sample results for PQL estimation in such a way so as to do not
depend on the value of $\hat{\bm{G}}$ itself [46, 40, in a spirit similar to
that of]; only the true random effects covariance matrix $\dot{\bm{G}}$
appears in our theorems.

We also adopt the above approach for the dispersion parameter in the GLMM.
That is, we set $\phi=\hat{\phi}$ in (2), where $\hat{\phi}$ is a known
constant or a stochastically bounded term. In the Poisson and binomial
distributions, $\hat{\phi}$ is set to its known true value $\dot{\phi}=1$. In
cases where the true dispersion parameter is unknown, we can use either a
constant or one of the suggested estimators of the dispersion parameter in the
literature (e.g., a scaled sum of squared conditional Pearson residuals). For
the remainder of this article, and as discussed in Section 1, we focus on the
fixed and random effects in GLMMs. We do not discuss inferential properties of
$\dot{\phi}$ and $\dot{\bm{G}}$.

##  3 Conditional on Random Effects

In many applications of independent-cluster GLMMs e.g., for longitudinal data,
covariates included as random effects are also included as fixed effects [38].
With this in mind, we develop our results under the setting where all
covariates are partnered i.e. included as both fixed and random effects such
that $\bm{x}_{ij}=\bm{z}_{ij}$ for all $(i,j)$ and $p_{f}=p_{r}=:p$. Next, let
$\bm{A}$ be a $q\times(m+1)p$ matrix with the finite selection property. That
is, for any row of $\bm{A}$, there exists an $m_{0}\in\mathbb{N}$ such that
the $\\{(m_{0}+1)p+1\\}$th to $\\{(m+1)p\\}$th components of the row are zero
for all $m>m_{0}$. All components of $\bm{A}$ must have a component-wise
limit, with at least one of these limits being non-zero. We partition $\bm{A}$
into $[\bm{A}_{f},\bm{A}_{r}]$, where $\bm{A}_{f}$ is of dimension $q\times p$
and $\bm{A}_{r}$ is of dimension $q\times mp$. Also, for an arbitrary matrix
$\bm{C}$, let $\bm{C}_{[i:j,k:l]}$ denote the sub-matrix comprising the $i$th
to $j$th row and $k$th to $l$th column of $\bm{C}$ and $\bm{C}_{[i,]}$ and
$\bm{C}_{[,j]}$ denote the $i$th row and $j$th column respectively. Similarly,
for a vector $\bm{c}$ we let $\bm{c}_{[i:j]}$ denote the sub-vector formed by
taking the $i$th to $j$th components; the quantity $\bm{c}_{[i]}$ simply
denotes the $i$th component of $\bm{c}$.

Let
$\bm{\mu}_{i}(\bm{\theta})=\\{a^{\prime}(\eta_{i1}),\ldots,a^{\prime}(\eta_{in_{i}})\\}^{\top}$,
$\bm{\mu}(\bm{\theta})=\\{a^{\prime}(\eta_{11}),\ldots,a^{\prime}(\eta_{mn_{m}})\\}^{\top}$,
$\dot{\bm{W}}_{i}=\hat{\phi}^{-1}\text{diag}\\{a^{\prime\prime}(\dot{\eta}_{i1}),\ldots,a^{\prime\prime}(\dot{\eta}_{in_{i}})\\}$
and
$\dot{\bm{W}}=\hat{\phi}^{-1}\text{diag}\\{a^{\prime\prime}(\dot{\eta}_{11}),\ldots,a^{\prime\prime}(\dot{\eta}_{mn_{m}})\\}$.
Furthermore, write $\dot{\mu}_{ij}=a^{\prime\prime}(\dot{\eta}_{ij})$,
$\dot{\bm{\mu}}_{i}=\bm{\mu}_{i}(\dot{\bm{\theta}})$ and
$\dot{\bm{\mu}}=\bm{\mu}(\dot{\bm{\theta}})$, and let $\otimes$ denote the
Kronecker product operator, $\bm{I}_{m}$ denote the $m\times m$ identity
matrix, and $\bm{1}_{m}$ denote a matrix or vector of ones, with dimension
indicated by the relevant subscripts. Furthermore, let
$\bm{D}_{r}=\text{diag}({n}_{1}^{1/2}\bm{1}_{p},\ldots,{n}_{m}^{1/2}\bm{1}_{p})$,
$\bm{D}=\text{bdiag}({N}^{1/2}\bm{I}_{p},\bm{D}_{r})$,
$\bm{D}^{*}=\text{bdiag}({m}^{1/2}\bm{I}_{p},\bm{D}_{r})$,
$\bm{D}^{+}={m}^{1/2}\bm{I}_{(m+1)p}$, and define the two limiting quantities

$\displaystyle\bm{\Omega}$
$\displaystyle=\underset{m,n_{L}\rightarrow\infty}{\lim}\frac{\dot{\phi}}{\hat{\phi}}\bm{A}\,\text{bdiag}\Bigg{\\{}\frac{1}{m}\sum_{i=1}^{m}\frac{n}{n_{i}}\left(\frac{\bm{X}^{\top}_{i}\dot{\bm{W}}_{i}\bm{X}_{i}}{n_{i}}\right)^{-1},\left(\frac{\bm{X}^{\top}_{1}\dot{\bm{W}}_{1}\bm{X}_{1}}{n_{1}}\right)^{-1},\ldots,\left(\frac{\bm{X}^{\top}_{m}\dot{\bm{W}}_{m}\bm{X}_{m}}{n_{m}}\right)^{-1}\Bigg{\\}}\bm{A}^{\top},$
$\displaystyle\bm{\Omega}_{r}$
$\displaystyle=\lim_{m,n_{L}\rightarrow\infty}\frac{\dot{\phi}}{\hat{\phi}}\bm{A}_{r}\bm{D}_{r}\left(\bm{Z}^{\top}\dot{\bm{W}}\bm{Z}\right)^{-1}\bm{D}_{r}^{\top}\bm{A}_{r}^{\top}.$

Note $\bm{\Omega}$ and $\bm{\Omega}_{r}$ are not actually functions of
$\hat{\phi}$, since
$\hat{\phi}\dot{\phi}^{-1}\dot{\bm{W}}_{i}=\dot{\phi}^{-1}\text{diag}\\{a^{\prime\prime}(\dot{\eta}_{i1}),\ldots,a^{\prime\prime}(\dot{\eta}_{in_{i}})\\}$
and similarly for $\dot{\bm{W}}$.

We consider the setting where both the minimum cluster size $n_{L}$ and the
number of clusters $m$ grow to infinity, such that $n_{i}=O(n_{L})$ uniformly
for $i=1,\ldots,m$. This implies for any $i=1,\ldots,m$, we have $n_{i}=O(n)$,
$n=O(n_{i})$, $N=O(mn_{i})$, and $mn_{i}=O(N)$. This restriction on the growth
rates of the cluster sizes is commonly employed in asymptotic analysis of PQL
estimation [60, e.g.,]. Additionally, we require the following regularity
conditions.

  1. (C1)

The function $a(\eta)$ is at least three times continuously differentiable in
its domain, with $0<c_{0}\leq a^{\prime\prime}(\eta)\leq c_{0}^{-1}<\infty$
and $|a^{\prime\prime\prime}(\eta)|\leq c_{0}^{-1}<\infty$ for some
sufficiently small constant $c_{0}$.

  2. (C2)

For every $i=1,\ldots,m$ and $j=1,\ldots,n_{i}$, there exists a sufficiently
large constant $C_{1}$ such that $\|\bm{x}_{ij}\|_{\infty}<C_{1}$ where
$\|\cdot\|_{\infty}$ is the maximum norm. Furthermore, denote
$\dot{\bm{H}}_{i}=(n_{i}^{-1}\hat{\phi}\dot{\phi}^{-1}\bm{X}_{i}^{\top}\dot{\bm{W}}_{i}\bm{X}_{i})^{-1}$.
Then for all $i=1,\ldots,m$, the matrices
$\lim_{n_{i}\rightarrow\infty}\dot{\bm{H}}_{i}=\dot{\bm{K}}_{i}$ and
$\lim_{m,n_{L}\rightarrow\infty}m^{-1}\sum_{i=1}^{m}nn_{i}^{-1}\dot{\bm{H}}_{i}=\dot{\bm{K}}$
are positive definite with minimum and maximum eigenvalues bounded from above
and below by $c_{1}^{-1}$ and $c_{1}$ respectively, for a sufficiently small
constant $c_{1}$.

  3. (C3)

The vector of true parameters
$\dot{\bm{\theta}}=(\dot{\bm{\beta}}^{\top},\dot{\bm{b}}^{\top})^{\top}$,
where
$\dot{\bm{b}}=(\dot{\bm{b}}_{1}^{\top},\ldots,\dot{\bm{b}}_{m}^{\top})^{\top}$,
is an interior point in some compact set $\Theta\subset\mathbb{R}^{(m+1)p}$.

  4. (C4)

The working matrix $\hat{\bm{G}}$ is positive definite, and its inverse is
$O_{p}(1)$. Also, the working quantity $\hat{\phi}$ is strictly positive and
$O_{p}(1)$.

  5. (C5)

For all $i=1,\ldots,m$ and $n_{i}\in\mathbb{N}$, it holds that
$E([n_{i}^{1/2}(\bm{X}_{i}^{\top}\dot{\bm{W}}_{i}\bm{X}_{i}+\hat{\bm{G}}^{-1})^{-1}\\{\hat{\phi}^{-1}\bm{X}_{i}^{\top}(\bm{y}_{i}-\dot{\bm{\mu}}_{i})-\hat{\bm{G}}^{-1}\dot{\bm{b}}_{i}\\}]^{4})<\infty$,
where the power and expectation are applied component-wise.

Conditions (C1) - (C3) are needed to guarantee the existence and regular
behavior of the asymptotic variance for the PQL estimating function, and to
establish a Lindeberg condition needed to obtain a central limit theorem.
Condition (C4) is required to ensure that the shrinkage of the random effects
is asymptotically negligible, and formalises our discussion of $\hat{\bm{G}}$
and $\hat{\phi}$ at the end of Section 2. Condition (C5) is needed to bound
the order of $\|\hat{\bm{\theta}}-\dot{\bm{\theta}}\|_{\infty}$, and is
satisfied by many distributions e.g. Poisson and binomial, when the random
effects are normally distributed [59, see also].

For the remainder of this section, we consider the regime where we condition
on the random effects, so that $\dot{\bm{\theta}}$ is a $(m+1)p$-vector of
constants. The assumptions and conditions outlined above however will be
applied to both the conditional and unconditional regime.

###  3.1 Main Result for the Conditional Regime

Let $\bm{1}^{*}_{m}=(-1,\bm{1}_{m}^{\top})^{\top}$. Then we have the
following:

######  Theorem 1.

Assume Conditions (C1) - (C5) are satisfied and $mn_{L}^{-1}\rightarrow 0$.
Then as $m,n_{L}\rightarrow\infty$, and conditional on the true vector of
random effects $\dot{\bm{b}}$, it holds that

  1. (a)

$\|\hat{\bm{\theta}}-(\dot{\bm{\theta}}-\bm{1}^{*}_{m}\otimes
m^{-1}\sum_{i=1}^{m}\dot{\bm{b}}_{i})\|_{\infty}=o_{p}(1)$.

  2. (b)

$\bm{A}\bm{D}\\{\hat{\bm{\theta}}-(\dot{\bm{\theta}}-\bm{1}^{*}_{m}\otimes
m^{-1}\sum_{i=1}^{m}\dot{\bm{b}}_{i})\\}\overset{D}{\rightarrow}N(\bm{0},\bm{\Omega})$.

The first part of the theorem establishes consistency for the PQL estimator
around a sum-to-zero reparametrized version of the true parameters (see below
for more discussion on the latter aspect). The block diagonal structure of
$\bm{\Omega}$ in the second part of the theorem shows that conditional on true
random effects vector, the corresponding estimators are asymptotically
independent between clusters, and also asymptotically independent of the fixed
effects estimators.

We illustrate a few special cases of Theorem 1 using specific selection
matrices. First, suppose $\bm{A}=[\bm{I}_{p},\bm{0}_{p\times mp}]$. If
$\sum_{i=1}^{m}\dot{\bm{b}}_{i}=\bm{0}_{p}$, then we obtain
$\bm{A}\bm{D}(\hat{\bm{\theta}}-\dot{\bm{\theta}})={N}^{1/2}(\hat{\bm{\beta}}-\dot{\bm{\beta}})\overset{D}{\rightarrow}N(\bm{0},\dot{\bm{K}})$
conditional on the random effects, where $\dot{\bm{K}}$ is the limiting matrix
defined in Condition (C2). Also, suppose
$\bm{A}=[\bm{0}_{p},\bm{I}_{p},\bm{0}_{p\times(m-1)p}]$. Then from Theorem 1,
we have
$\bm{A}\bm{D}(\hat{\bm{\theta}}-\dot{\bm{\theta}})={n}_{1}^{1/2}(\hat{\bm{b}}_{1}-\dot{\bm{b}}_{1})\overset{D}{\rightarrow}N(\bm{0},\dot{\bm{K}}_{1})$,
conditional on the random effects. The analogous result holds for choosing any
particular cluster. Finally, since the random effects exhibit a slower
convergence rate than the fixed effects, and noting the asymptotic
independence, then for an arbitrary $p$-dimensional constant $\bm{a}$ we
obtain
${n}_{i}^{1/2}\bm{a}^{\top}(\hat{\bm{\beta}}+\hat{\bm{b}}_{i}-\dot{\bm{\beta}}-\dot{\bm{b}}_{i})\overset{D}{\rightarrow}N\left(\bm{0},\bm{a}^{\top}\dot{\bm{K}}_{i}\bm{a}\right);\quad
i=1,\ldots,m$, conditional on the random effects. As an example, if the linear
predictor involves a fixed and random intercept and a fixed and random slope
for a single covariate, then we set $\bm{a}=(1,x_{ij})^{\top}$ and obtain
$n_{i}^{1/2}(\hat{\eta}_{ij}-\dot{\eta}_{ij})=n_{i}^{1/2}(\hat{\beta}_{0}+\hat{b}_{i0}+\hat{\beta}_{1}x_{ij}+\hat{b}_{i1}x_{ij}-\dot{\beta}_{0}-\dot{b}_{i0}-\dot{\beta}_{1}x_{ij}-\dot{b}_{i1}x_{ij})\overset{d}{\rightarrow}N(0,\bm{a}^{\top}\dot{\bm{K}}_{i}\bm{a})$.

For statistical inference, we can appeal to Slutsky’s Theorem and replace
$\dot{\bm{K}}_{i}$ with $\hat{\bm{H}}_{i}$, and $\dot{\bm{K}}$ with
$m^{-1}\sum_{i=1}^{m}nn_{i}^{-1}\hat{\bm{H}}_{i}$. Here $\hat{\bm{H}}_{i}$ is
defined as $(n_{i}^{-1}\bm{X}_{i}^{\top}\hat{\bm{W}}_{i}\bm{X}_{i})^{-1}$
where
$\hat{\bm{W}}_{i}=\tilde{\phi}^{-1}\text{diag}\\{a^{\prime\prime}(\hat{\eta}_{i1}),\ldots,a^{\prime\prime}(\hat{\eta}_{in_{i}})\\}$
for some consistent estimator of the dispersion parameter $\tilde{\phi}$ e.g.,
based on the inverse scaled sum of squared conditional Pearson residuals.
Theorem 1 then provides a straightforward way to construct confidence
intervals, say, for all the parameters and combinations thereof. In fact, the
forms of these intervals are similar to standard results in (penalized)
GLMs[52]: this is not surprising given we are working conditional on the true
vector of random effects. say.

Finally, note the PQL estimator is consistent for a sum-to-zero reparametrized
version of the true parameters. This occurs because the PQL estimators of the
random effects must satisfy a sum-to-zero constraint regardless of the
underlying true parameter values, and under a conditional regime, this induces
an asymptotic bias
$\bm{1}^{*}_{m}\otimes(m^{-1}\sum_{i=1}^{m}\dot{\bm{b}}_{i})$ in Theorem 1,
which can be interpreted as shifting the mean of the random effects into the
corresponding fixed effects. We offer more discussion around this asymptotic
bias in the supplementary material.

##  4 Unconditional Regime

We now turn to establishing results under an unconditional regime i.e.,
treating $\dot{\bm{b}}_{i}$’s as random instead of conditioning on them. This
has two main implications. First, in the unconditional setting the quantity
$m^{-1}\sum_{i=1}^{m}\dot{\bm{b}}_{i}$ is no longer deterministic and should
not be treated as a bias term. Instead, it is of order $O_{p}(m^{-1/2})$, and
so competes with other leading terms in the relevant Taylor expansion to be
the dominating term. This results in a reduction of the rate of convergence
for the fixed effects estimator, from $N^{1/2}$ in the conditional regime to
$m^{1/2}$ in the unconditional regime. Second, in contrast to the conditional
regime, the observations within the same cluster are no longer independent.
This has ramifications when applying the central limit theorem to establish
asymptotic multivariate normality. In Section 4.1, we provide a simple but
insightful example based on a Poisson random intercept model, which
demonstrates that the prediction gap is not always asymptotically normally
distributed.

The two approximations below, derived from the Taylor expansion of the PQL
objective function, will be central to understanding the large sample
developments we make on a more intuitive level. For a given $\hat{\phi}$, we
have

$\displaystyle\hat{\bm{\beta}}-\dot{\bm{\beta}}$
$\displaystyle=m^{-1}\sum_{i=1}^{m}\dot{\bm{b}}_{i}+o_{p}(1)$ (3a)
$\displaystyle\hat{\bm{b}}-\dot{\bm{b}}$ $\displaystyle=-\bm{1}_{m}\otimes
m^{-1}\sum_{i=1}^{m}\dot{\bm{b}}_{i}+(\bm{Z}^{\top}\dot{\bm{W}}\bm{Z})^{-1}\\{{\hat{\phi}}^{-1}\bm{Z}^{\top}(\bm{y}-\dot{\bm{\mu}})\\}+o_{p}(1).$
(3b)

We will refer to both equations in the discussion of the theorems to be
presented later on.

###  4.1 Prediction Gap - Counterexample

We offer a motivating and insightful example to illustrate that the prediction
gap is not, in general, asymptotically normally distributed. This example also
offers a simple case where $\bm{X}_{i}\neq\bm{Z}_{i}$, and offers an
interesting comparison to the theory established under the assumption of
$\bm{X}_{i}=\bm{Z}_{i}$.

Consider a Poisson random intercept model with canonical log link. That is,
the true model is given by
$f(y_{ij}|\dot{b}_{i})=\exp(y_{ij}\dot{\eta}_{ij}-\dot{\mu}_{ij})/(y_{ij}!)$
with $\ln(\dot{\mu}_{ij})=\dot{\eta}_{ij}=\dot{b}_{i}$, and
$\dot{b}_{i}\overset{i.i.d.}{\sim}N(0,\dot{\sigma}_{b}^{2})$. Assume a working
$\hat{\sigma}_{b}^{2}$, and apply PQL estimation to estimate the random
effects $b_{i}$ for $i=1,\ldots,n$. For simplicity, we also assume a balanced
design, such that $n_{i}=n$ for all $i=1,\ldots,m$. Then it is possible to
show (see the supplementary material for the formal derivation) that when
$mn^{-2}\rightarrow 0$, the prediction gap of the first cluster
$\hat{b}_{1}-\dot{b}_{1}$ satisfies

$\displaystyle{n}^{1/2}(\hat{b}_{1}-\dot{b}_{1})$
$\displaystyle=n^{-1/2}\sum_{j=1}^{n}\\{y_{1j}\exp(-\dot{b}_{1})-1\\}+o_{p}(1).$
(4)

Therefore, we obtain $\hat{b}_{1}=\dot{b}_{1}+o_{p}(1)$, and similarly for
each cluster $i=1,\ldots,m$. When conditioned on $\dot{b}_{1}$,
$n^{-1/2}\sum_{j=1}^{n}\\{y_{1j}\exp(-\dot{b}_{1})-1\\}$ is a normalised sum
of independent random variables. Unconditionally however, the sum consists of
an exchangeable collection of uncorrelated but dependent random variables with
mean zero and finite non-zero variance. Using the central limit theorem for
exchangeable random variables [33], it can be subsequently be shown that the
quantity $n^{-1/2}\sum_{j=1}^{n}\\{y_{1j}\exp(-\dot{b}_{1})-1\\}$, and hence
${n}^{1/2}(\hat{b}_{1}-\dot{b}_{1})$, is not asymptotically normally
distributed.

With the above example in mind, we now state the main results for the
unconditional regime.

###  4.2 Fixed Effects

We have the following result for the PQL estimator of the fixed effects under
an unconditional regime.

######  Theorem 2.

Assume Conditions (C1) - (C5) are satisfied, and $mn_{L}^{-2}\rightarrow 0$.
Then as $m,n_{L}\rightarrow\infty$ and unconditional on the random effects
$\dot{\bm{b}}$, it holds that
${m}^{1/2}(\hat{\bm{\beta}}-\dot{\bm{\beta}})\overset{D}{\rightarrow}N(\bm{0},\dot{\bm{G}})$.

This result should not be too surprising given the form of (3a). Furthermore,
the rate of convergence and asymptotic distribution coincides with the result
obtained by [47] for the partnered fixed effects for the quasi-maximum
likelihood estimator. More importantly, Theorem 2 allows practitioners to
straightforwardly perform statistical inference for the fixed effects, so long
as $mn_{L}^{-2}\rightarrow 0$. Although $\dot{\bm{G}}$ is not known, we can
appeal to Slutsky’s theorem and replace it with a consistent estimator (e.g.,
the sample covariance matrix of the estimated random effects). Theorem 2
contrasts with Theorem 1 derived under the conditional regime, where
$mn_{L}^{-1}\rightarrow 0$ is required but the convergence rate is
${N}^{1/2}$. This reduction in the rate of convergence arises because the
leading term in the Taylor expansion is different: in the unconditional
regime, it is simply the normalised sum of random effects over all the
clusters, and thus its variability is dominated by the term
$m^{-1/2}\sum_{i=1}^{m}\dot{\bm{b}}_{i}$. However, this term is deterministic
in the conditional regime, and serves to enforce a sum-to-zero constraint
instead as discussed in Section 3. Generally speaking, the Taylor expansion
can be interpreted as comprising terms which either capture the stochasticity
in the random effects vector $\dot{\bm{b}}$, or the stochasticity in responses
$y_{ij}$ given the random effects. These terms compete with each other, and
which one dominates depends on the relative rates of $m$ and $n_{i}$. This
intricacy in the nature of the results will be made apparent in our results
for the prediction gap in Section 4.4.

###  4.3 Estimators of the Random Effects

Next, we state a convergence result for the PQL estimators of the random
effects under the unconditional regime.

######  Theorem 3.

Assume Conditions (C1) - (C5) are satisfied and $mn_{L}^{-2}\rightarrow 0$.
Then as $m,n_{L}\rightarrow\infty$ and unconditional on the random effects
$\dot{\bm{b}}$, it holds that
$\bm{A}_{r}(\hat{\bm{b}}-\dot{\bm{b}})\overset{P}{\rightarrow}\bm{0}_{q}$ .

Practically, Theorem 3 confirms the asymptotic distribution of a finite subset
of the PQL estimators is the distribution of the random effects themselves.
This can play a useful role for helping to validate the examination of the
empirical distribution of PQL estimators $\hat{\bm{b}}$ as a model diagnostic
tool. For instance, if the random effects are normally distributed and
$\bm{A}_{r}$ only selects the first cluster, then we would expect
$\hat{\bm{b}}_{1}$ to have an approximate $N(\bm{0},\dot{\bm{G}})$
distribution. On the other hand, Theorem 3 does not help us when it comes to
performing likelihood-based inference for the true random effects
$\dot{\bm{b}}$, as this does not appear in the approximation
$\hat{\bm{b}}_{1}\sim N(\bm{0},\dot{\bm{G}})$ itself.

As an aside, note the above means we can apply the continuous mapping theorem
and show that
$q^{-1}\sum_{i=1}^{q}\hat{\bm{b}}_{i}\hat{\bm{b}}_{i}^{\top}-q^{-1}\sum_{i=1}^{q}\dot{\bm{b}}_{i}\dot{\bm{b}}_{i}^{\top}\overset{P}{\rightarrow}0$
for any $q\in\mathbb{N}$. Since
$q^{-1}\sum_{i=1}^{q}\dot{\bm{b}}_{i}\dot{\bm{b}}_{i}^{\top}\overset{P}{\rightarrow}\dot{\bm{G}}$
as $q\rightarrow\infty$, this further reiterates the use of a sample
covariance matrix of the estimated random effects as an estimator of
$\dot{\bm{G}}$ [46, 43, consistent with].

###  4.4 Prediction Gap

In this section, we present a result for the large sample distribution of a
finite subset of the prediction gap, $\hat{\bm{b}}-\dot{\bm{b}}$, in the
unconditional regime. As mentioned above, the asymptotic distribution as well
as the convergence rate of the prediction gap depends on the relative rates of
growth of $m$ and $n_{i}$. This contrasts with the conditional regime, where
there is no dependence on the relative rate and the PQL estimator of the
random effects is always normally distributed with the convergence rate
$n_{i}^{1/2}$.

We first introduce some terminology. Suppose we have two arbitrary continuous
cumulative distribution functions (cdfs) $F_{1}$ and $F_{2}$ with supports in
$\mathbb{R}^{p}$. Then we define the convolution of $F_{1}$ and $F_{2}$,
denoted $F_{1}*F_{2}$, as
$(F_{1}*F_{2})(\bm{z})=\int_{\mathbb{R}^{p}}F_{1}(\bm{z}-\bm{\tau})dF_{2}(\bm{\tau})$.
Next, for a random variable $\bm{P}$, we say
$\bm{P}\sim\text{mixN}\\{\bm{\mu}(\bm{b}),\bm{\Sigma}(\bm{b}),F_{\bm{b}}\\}$
if $\bm{P}|\bm{b}\sim N\\{\bm{\mu}(\bm{b}),\bm{\Sigma}(\bm{b})\\}$ and
$F_{\bm{b}}$ is the cdf of $\bm{b}$, where the conditional mean vector
$\bm{\mu}(\bm{b})$ and covariance matrix $\bm{\Sigma}(\bm{b})$ may depend on
$\bm{b}$. In other words, $\bm{P}$ has cdf
$F_{\bm{P}}(\bm{p})=\int\Psi_{\bm{P}|\bm{b}}(\bm{p})dF_{\bm{b}}(\bm{b})$,
where $\Psi_{\bm{P}|\bm{b}}$ is the cdf of
$N\\{\bm{\mu}(\bm{b}),\bm{\Sigma}(\bm{b})\\}$. A special case of this normal
scale-mixture distribution is when $\bm{\mu}(\bm{b})$ and
$\bm{\Sigma}(\bm{b})$ do not depend on $\bm{b}$, so that
$F_{\bm{P}}(\bm{p})=\int\Psi_{\bm{P}|\bm{b}}(\bm{p})dF_{\bm{b}}(\bm{b})=\Psi_{\bm{P}|\bm{b}}(\bm{p})\int
dF_{\bm{b}}(\bm{b})=\Psi_{\bm{P}|\bm{b}}(\bm{p})$; in other words, the normal
scale-mixture distribution reduces to a normal distribution. Note estimators
with asymptotic normal mixture distributions have arisen in previous
literature, for instance, on results relating to local asymptotic normality
and non-ergodic models [37, 31].

Using the above definition, we obtain the following results.

######  Theorem 4.

Assume Conditions (C1)-(C5) are satisfied and $mn_{L}^{-2}\rightarrow 0$. Then
as $m,n_{L}\rightarrow\infty$ and unconditional on the random effects
$\dot{\bm{b}}$, for each $i=1,\ldots,m$ we have the following:

  1. (a)

If $mn_{i}^{-1}\rightarrow\infty$, then
${n}_{i}^{1/2}(\hat{\bm{b}}_{i}-\dot{\bm{b}}_{i})\overset{D}{\rightarrow}\text{mixN}(\bm{0},\dot{\bm{K}}_{i},F_{\dot{\bm{b}}_{i}})$.

  2. (b)

If $mn_{i}^{-1}\rightarrow\gamma_{i}\in(0,\infty)$, then
${n}_{i}^{1/2}(\hat{\bm{b}}_{i}-\dot{\bm{b}}_{i})\overset{D}{\rightarrow}\text{mixN}(\bm{0},\dot{\bm{K}}_{i},F_{\dot{\bm{b}}_{i}})*N(\bm{0},\gamma_{i}^{-1}\dot{\bm{G}})$.

  3. (c)

If $mn_{i}^{-1}\rightarrow 0$, then
${m}^{1/2}(\hat{\bm{b}}_{i}-\dot{\bm{b}}_{i})\overset{D}{\rightarrow}N(\bm{0},\dot{\bm{G}})$.

######  Corollary 1.

Assume Conditions (C1)-(C5) are satisfied, and $mn_{L}^{-2}\rightarrow 0$. If
$mn_{L}^{-1}\rightarrow\infty$, then as $m,n_{L}\rightarrow\infty$ and
unconditional on the random effects $\dot{\bm{b}}$,
$\bm{A}_{r}\bm{D}_{r}(\hat{\bm{b}}-\dot{\bm{b}})\overset{D}{\rightarrow}\text{mixN}(\bm{0},\bm{\Omega}_{r},F_{\dot{\bm{b}}})$.

Theorems 3 and 4 bears some similarity to the results of [50], who show for
LMMs that the distribution of the EBLUP can asymptotically be written as the
convolution between the distribution of the random effects and the
distribution of a smaller order stochastic term. However, the above is the
first to establish such results for GLMMs. Theorem 4 states that the correct
asymptotic distribution to use when performing inference using the PQL
estimate of the random effects depends on the relative growth rates of $m$ and
$n_{i}$. As hinted at previously, this is a consequence of there being two
competing terms in the corresponding Taylor expansion (3b): one term arising
from the random effects, and the other term arising from the distribution of
the responses given the random effects.

When $mn_{i}^{-1}\rightarrow\infty$ i.e., the number of clusters grows faster
than the cluster size, the appropriate asymptotic distribution is given by the
scale-mixture distribution
$\text{mixN}\\{\bm{0},(\hat{\phi}\dot{\phi}^{-1}\bm{X}_{i}^{\top}\dot{\bm{W}}_{i}\bm{X}_{i})^{-1},F_{\dot{\bm{b}}_{i}}\\}$,
noting again that
$\hat{\phi}\dot{\phi}^{-1}\dot{\bm{W}}_{i}=\dot{\phi}^{-1}\text{diag}\\{a^{\prime\prime}(\dot{\eta}_{i1}),\ldots,a^{\prime\prime}(\dot{\eta}_{in_{i}})\\}$.
Corollary 1 offers a slightly more general result than that given in Theorem 4
for the $mn_{L}^{-1}\rightarrow\infty$ case. Note in the linear case, the GLM
iterative weights $\dot{\bm{W}}$ do not depend on the random effects
$\dot{\bm{b}}$, and so the corresponding normal scale-mixture distribution
reduces to a normal distribution, consistent with the asymptotic normality
result derived for the EBLUP in [50]. Practically, numerical techniques or
simulation are required to compute the quantiles of the normal scale-mixture
distribution for constructing prediction intervals. We use this approach in
our simulations in Section 5.

When $mn_{i}^{-1}\rightarrow 0$ i.e., the cluster sizes grows faster than the
number of clusters, Theorem 4 shows that the appropriate approximation to
consider is the normal distribution $N(\bm{0},n^{-1}\dot{\bm{G}})$. Note this
is identical to the fixed effects result of Theorem 2, and yields relatively
straightforward prediction intervals for $\dot{\bm{b}}_{i}$ as long as we have
a consistent estimator for $\dot{\bm{G}}$. Intuitively, the asymptotic
distribution here is identical to that derived in Theorem 2 because the
dominating terms in the Taylor expansions in both cases are effectively the
same. Finally, when $mn_{i}^{-1}\rightarrow\gamma\in(0,\infty)$, Theorem 4b
states that the asymptotic distribution of the PQL estimates is given by the
convolution of the two cases above, noting that these two leading terms in the
Taylor expansion are asymptotically independent. Again, numerical
techniques/simulations are needed to compute prediction intervals.

In summary, Theorem 4 offers an asymptotically valid way of computing
prediction intervals for the realised random effects in the unconditional
regime, when the random effects have a corresponding partnered fixed effect in
the model. It implies that estimating the variance of the prediction gap, and
then naively assuming normality in order to construct prediction intervals for
the random effects, will fail to yield asymptotically correct inference under
the unconditional regime for PQL estimation.

###  4.5 Linear Predictor

Neither Theorems 2 nor 4 above derive the joint distribution of the fixed
effects estimator and prediction gap, of which the linear predictor is a
function. Below, to address this, we establish a separate result specifically
for the sum of a random effect and its partnered fixed effect, given an
arbitrary $p$-dimensional constant vector $\bm{a}$.

######  Theorem 5.

Assume Conditions (C1)-(C5) are satisfied, $mn_{L}^{-2}\rightarrow 0$, and
$mn_{U}^{-1/2}\rightarrow\infty$. Then as $m,n_{L}\rightarrow\infty$ and
unconditional on the random effects $\dot{\bm{b}}$, it holds for each
$i=1,\ldots,m$ that
${n}_{i}^{1/2}\bm{a}^{\top}(\hat{\bm{\beta}}+\hat{\bm{b}}_{i}-\dot{\bm{\beta}}-\dot{\bm{b}}_{i})\overset{D}{\rightarrow}\text{mixN}(\bm{0},\bm{a}^{\top}\dot{\bm{K}}_{i}\bm{a},F_{\dot{\bm{b}}_{i}})$.

As an example, consider again a linear predictor involving a fixed and random
intercept and a fixed and random slope for a single covariate. Then we set
$\bm{a}=(1,x_{ij})^{\top}$ and obtain
$n_{i}^{1/2}(\hat{\eta}_{ij}-\dot{\eta}_{ij})=n_{i}^{1/2}(\hat{\beta}_{0}+\hat{b}_{i0}+\hat{\beta}_{1}x_{ij}+\hat{b}_{i1}x_{ij}-\dot{\beta}_{0}-\dot{b}_{i0}-\dot{\beta}_{1}x_{ij}-\dot{b}_{i1}x_{ij})\overset{d}{\rightarrow}mixN(0,\bm{a}^{\top}\dot{\bm{K}}_{i}\bm{a},F_{\dot{\bm{b}}_{i}})$.
For performing inference on the linear predictor in a GLMM or functions
thereof, Theorem 5 states that we again need to employ the normal scale-
mixture distribution. This result differs from the asymptotic normality of the
linear predictor derived under the conditional regime in Section 3. Note we
can also develop a similar result for the difference between the prediction
gaps of two clusters, and we refer to the supplementary material for details
of this result.

To conclude the section, we remark that Theorems 2-5 do not offer results on
the joint distribution of the prediction gap and fixed effects. However, we
know from the associated proof that the prediction gaps for each cluster are
asymptotically independent from each other as well as from the fixed effects
estimator when $mn_{U}^{-1}\rightarrow\infty$ and $mn_{L}^{-2}\rightarrow 0$,
and so a joint distribution result can be derived from this.

##  5 Simulation Study

We performed a numerical study to assess the usefulness of our asymptotic
results in finite samples. We simulated data from an independent-cluster GLMM
with five fixed and random effect covariates, considering Poisson and
Bernoulli responses, as follows. First, we set the first component of
$\bm{x}_{ij}=\bm{z}_{ij}$ equal to one to represent a fixed/random intercept.
The second and third components are simulated from a bivariate normal
distribution with mean zero and standard deviation one, with correlation equal
to 0.5. The fourth component is generated independently from a standard normal
distribution, and the last component is simulated from a Bernoulli
distribution with a probability of success equal to 0.5. Next, we set the
5-vector of true fixed effect coefficients to either
$\dot{\bm{\beta}}=(2,0.1,-0.1,0.1,0.1)^{\top}$ for Poisson responses, or
$\dot{\bm{\beta}}=(-0.1,0.1,-0.1,0.1,0.1)^{\top}$ for Bernoulli responses and
the $5\times 5$ random effects covariance matrix in both cases to
$\dot{\bm{G}}=\bm{I}_{5}$. Based on these true parameter values, we simulated
the random effect coefficients
$\dot{\bm{b}}_{i}\overset{i.i.d.}{\sim}\mathcal{N}(0,\dot{\bm{G}})$. Finally,
conditional on $\dot{\bm{b}}_{i}$ the responses $y_{ij}$ were generated from
either a Poisson distribution with log link, or a Bernoulli distribution with
logit link. We varied the number of clusters as $m=\\{25,50,100,200,400\\}$
and the cluster sizes $n_{i}=n=\\{25,50,100,200,400\\}$, noting we assumed
equal cluster sizes in the simulation design for simplicity For each
combination of $(m,n)$, we simulated 1000 datasets. For the conditional
regime, we simulated $\dot{\bm{b}}$ only once and conditioned on this for all
simulated datasets; for unconditional regime, we simulated a new
$\dot{\bm{b}}$ for each simulated dataset.

For each simulated dataset, we fitted the corresponding GLMM using PQL
estimation, where we use the sample covariance matrix of the estimated random
effects as our update for $\hat{\bm{G}}$. That is, we iteratively maximize
equation (2) with respect to $\bm{\beta}$ and $\bm{b}$ for a given
$\hat{\bm{G}}$ (noting $\hat{\phi}=1$ is known for both these distributions),
and update $\hat{\bm{G}}$ as
$m^{-1}\sum_{i=1}^{m}\hat{\bm{b}}_{i}\hat{\bm{b}}^{\top}_{i}$, until
convergence.

We assessed performance separately under the conditional and unconditional
regimes. In the former, we examined the empirical coverage probability of 95%
coverage intervals constructed for $\bm{\beta}$ and for $\bm{b}_{1}$ (the
choice of the first cluster is arbitrary). The intervals were constructed
based on Theorem 1, with the asymptotic covariance matrix $\bm{\Omega}$
computed using the true parameter values. We refer to such intervals as
coverage intervals as opposed to confidence intervals. We also performed
Shapiro-Wilk tests on the components of the (1000) realised PQL estimates of
$\bm{\beta}$ and $\bm{b}_{1}$, in order to assess the asymptotic normality of
their respective sampling distributions. For the unconditional regime, we
examined the empirical coverage probability of 95% coverage intervals
constructed from Theorems 2 and 4 respectively. Again, this was done for the
fixed effect coefficients $\bm{\beta}$ and the random effects for the the
first cluster $\bm{b}_{1}$. To construct all intervals, we used the true
parameter values to compute the relevant asymptotic variance (this was done
solely to reduce the computational burden of the numerical study), and, when
required, obtained quantiles of relevant normal scale-mixture distributions by
directly simulating 10,000 samples from them. We also performed Shapiro-Wilk
tests on the components of the (1000) realised values of
$\hat{\bm{\beta}}-\dot{\bm{\beta}}$ and $\hat{\bm{b}}_{1}-\dot{\bm{b}}_{1}$.
Finally, we examined histograms for the third components of
$\hat{\bm{\beta}}-\dot{\bm{\beta}}$ and $\hat{\bm{b}}_{1}-\dot{\bm{b}}_{1}$,
which are representative of the histograms of the other components, as an
additional method of assessing asymptotic normality of the corresponding
sampling distributions.

###  5.1 Simulation Results

For reasons of brevity, below we focus on results for the unconditional
regime. Results for the conditional regime are presented in the supplementary
material and largely support the use of Theorem 1 for inference.

For the unconditional regime, Figures 1 and 2 display the empirical coverage
probabilities and results from applying the Shapiro-Wilk test, respectively.
For the fixed effect coefficients, the coverage probabilities for the
intervals obtained based on Theorem 2 were relatively accurate across most
combinations of $(m,n)$, with the exception of when $(m,n)=(25,25)$. For the
random effect coefficients, the coverage probabilities for intervals
calculated based on Theorem 4 approached the nominal coverage rapidly as
$(m,n)$ increased for the Poisson response case, while for the Bernoulli case
convergence was slightly slower due to the reduced amount of information per
response.

Figure 1: Empirical coverage probability of 95% coverage intervals for the
fixed and random effects, obtained under the unconditional regime.

Figure 2: $p$-values from Shapiro-Wilk tests applied to the fixed and random
effects estimates obtained using maximum PQL estimation, under the
unconditional regime.

Figure 3: Histograms for the third components of
$\hat{\bm{\beta}}-\dot{\bm{\beta}}$ (left panels) and
$\hat{\bm{b}}_{1}-\dot{\bm{b}}_{1}$ (right panels), under the unconditional
regime. Vertical facets represent the cluster sizes, while horizontal facets
represent the number of clusters. The dotted blue line indicates zero, and the
red curve is a kernel density smoother.

The Shapiro-Wilk tests run were consistent with the conclusions reached in
Theorems 2–4. Specifically, PQL estimates of the fixed effect coefficients
generally did not exhibit signs of non-normality, but the _difference_ between
the estimators and true random effects displayed evidence of non-normality
except when $n$ grew faster than $m$. This is also supported by the histograms
in Figure 3 which show some evidence of higher kurtosis in the cases
corresponding to small p-values in the Shapiro Wilk test. The histograms also
suggest that both $m$ and $n$ need to grow for the estimators to be consistent
for the true fixed and random effects, and in particular $n$ needs to grow for
the estimators to be unbiased. This is true especially for the Bernoulli
responses, for which convergence was much slower and very large cluster sizes
were needed for the estimators to be relatively unbiased.

In the supplementary material, we present additional results which showed that
the sample covariance matrix of the estimated random effects became a better
estimator of the true random effects covariance matrix $\dot{\bm{G}}$ as both
$m$ and $n$ grew. Also, recall from our discussion in Section 2 that our
asymptotic developments only require a working $\hat{\bm{G}}$, which need not
be a consistent estimator of the true random effects covariance matrix. As a
demonstration of this, we performed several additional simulations where in
the PQL estimation procedure, we simply fix $\hat{\bm{G}}$ to a constant
matrix and considered choices e.g., some constant multiplied by the identity
matrix. Results in the supplementary material demonstrate that coverage
probabilities for our proposed intervals still tended to the nominal level as
$(m,n)$ increased, while corresponding Shapiro-Wilk tests and histograms were
also consistent with our theory in large sample sizes and the empirical
results presented above.

##  6 Discussion

In this article, we established new asymptotic distributional results for
fixed effects, random effects, and the prediction gap, for an independent-
cluster GLMM fitted using penalized quasi-likelihood estimation. Our results
have important implications when it comes to inference and prediction for
mixed-effects models. For the conditional regime, we establish asymptotic
normality for any finite subset of the parameters. For random effects
predictions and inference in the unconditional regime, we validate examining
the empirical distribution of the estimated random effects as a diagnostic
tool for assessing deviations away from the assumed random effects
distribution [44, as is already commonly done in practice for GLMMs e.g.,]. On
the other hand, while the random effects estimators obtained using PQL are
asymptotically normally distributed when the true random effects are normally
distributed, we demonstrate that the difference between these two i.e., the
prediction gap, need not be normally distributed. Our large sample results
thus suggest the use of a normal approximation when performing unconditional
inference for the random effects, as is commonly done in practice [32, 36],
can be potentially misleading.

An important avenue of future research is to establish rates of convergence,
especially in the unconditional regime, when $\bm{x}_{ij}$ contains both
$\bm{z}_{ij}$ plus additional components which are only included as purely
fixed effects in the model. In the supplementary material, we develop some
further results for such unpartnered fixed effects in the special cases of
LMMs and GLMs. In both these cases, we see the convergence rate improves from
$O_{p}(m^{1/2})$ to $O_{p}(N^{1/2})$, compared to the partnered fixed effects.
On the other hand, for random effects without a partnered fixed effect, it is
likely that the correct asymptotic distribution for the prediction gap will be
the normal scale-mixture irrespective of the relative rates of $m$ and
$n_{i}$, as we saw in the motivating counterexample. Also, relaxing the
canonical link assumption is an interesting and important extension to our
work; we conjecture that non-canonical links could be accounted for by
generalising the form of the GLM iterative weights, as is done in GLMs.

## References

  * [1] Ishwar V Basawa and David John Scott  “Asymptotic optimal inference for non-ergodic models”  Springer Science & Business Media, 2012
  * [2] Douglas Bates, Martin Mächler, Ben Bolker and Steve Walker  “Fitting Linear Mixed-Effects Models Using lme4”  In _Journal of Statistical Software_ 67.1, 2015, pp. 1–48  DOI: 10.18637/jss.v067.i01
  * [3] JR Blum, H Chernoff, M Rosenblatt and H Teicher  “Central limit theorems for interchangeable processes”  In _Canadian Journal of Mathematics_ 10 Cambridge University Press, 1958, pp. 222–229
  * [4] Norman E Breslow and David G Clayton  “Approximate inference in generalized linear mixed models”  In _Journal of the American statistical Association_ 88, 1993, pp. 9–25
  * [5] Norman E Breslow and Xihong Lin  “Bias correction in generalised linear mixed models with a single component of dispersion”  In _Biometrika_ 82, 1995, pp. 81–91
  * [6] Mollie E. Brooks et al.  “glmmTMB Balances Speed and Flexibility Among Packages for Zero-inflated Generalized Linear Mixed Modeling”  In _The R Journal_ 9, 2017, pp. 378–400
  * [7] Lucien Le Cam and Grace L. Yang  “On the Preservation of Local Asymptotic Normality under Information Loss”  In _The Annals of Statistics_ 16 Institute of Mathematical Statistics, 1988, pp. 483–520
  * [8] Jing Cheng et al.  “Real longitudinal data analysis for real people: building a good enough mixed model”  In _Statistics in medicine_ 29, 2010, pp. 504–520
  * [9] Peter J Downey  “Distribution-free bounds on the expectation of the maximum with scheduling applications”  In _Operations Research Letters_ 9, 1990, pp. 189–201
  * [10] Yingying Fan and Runze Li  “Variable selection in linear mixed effects models”  In _Annals of statistics_ 40, 2012, pp. 2043–2068
  * [11] Charles R Henderson  “Sire evaluation and genetic trends”  In _Journal of Animal Science_ 1973, 1973, pp. 10–41
  * [12] Francis K C Hui  “On the use of a penalized quasilikelihood information criterion for generalized linear mixed models”  In _Biometrika_ 108, 2020, pp. 353–365
  * [13] Francis K C Hui, Samuel Müller and A H Welsh  “Joint selection in mixed models using regularized PQL”  In _Journal of the American Statistical Association_ 112, 2017, pp. 1323–1333
  * [14] Francis K C Hui, Samuel Müller and Alan H Welsh  “Random effects misspecification can have severe consequences for random effects inference in linear mixed models”  In _International Statistical Review_ 89, 2021, pp. 186–206
  * [15] Jiming Jiang  “Empirical best prediction for small-area inference based on generalized linear mixed models”  In _Journal of Statistical Planning and Inference_ 111.1-2 Elsevier, 2003, pp. 117–127
  * [16] Jiming Jiang, Haomiao Jia and Hegang Chen  “MAXIMUM POSTERIOR ESTIMATION OF RANDOM EFFECTS IN GENERALIZED LINEAR MIXED MODELS”  In _Statistica Sinica_ 11.1 Institute of Statistical Science, Academia Sinica, 2001, pp. 97–120
  * [17] Jiming Jiang, Matt P. Wand and Aishwarya Bhaskaran  “Usable and precise asymptotics for generalized linear mixed model analysis and design”  In _Journal of the Royal Statistical Society: Series B_ 84, 2022, pp. 55–82
  * [18] Raghu N Kackar and David A Harville  “Approximations for standard errors of estimators of fixed and random effects in mixed linear models”  In _Journal of the American Statistical Association_ 79 Taylor & Francis, 1984, pp. 853–862
  * [19] Łukasz Kidziński, Francis KC Hui, David I Warton and Trevor J Hastie  “Generalized Matrix Factorization: efficient algorithms for fitting generalized linear latent variable models to large data arrays”  In _The Journal of Machine Learning Research_ 23, 2022, pp. 13211–13239
  * [20] Ziyang Lyu and A H Welsh  “Asymptotics for EBLUPs: Nested Error Regression Models”  In _Journal of the American Statistical Association_ 117, 2021, pp. 1–15
  * [21] Ziyang Lyu and A H Welsh  “Increasing Cluster Size Asymptotics for Nested Error Regression Models”  In _Journal of Statistical Planning and Inference_ 217, 2021, pp. 52–68
  * [22] Charles E McCulloch and Shayle R Searle  “Generalized, linear, and mixed models”  John Wiley & Sons, 2004
  * [23] Lei Nie  “Convergence rate of MLE in generalized linear and nonlinear mixed-effects models: Theory and applications”  In _Journal of Statistical Planning and Inference_ 137, 2007, pp. 1787–1804
  * [24] HE Ogden  “On asymptotic validity of naive inference with an approximate likelihood”  In _Biometrika_ 104, 2017, pp. 153–164
  * [25] Helen Ogden  “On the error in Laplace approximations of high-dimensional integrals”  In _Stat_ 10, 2021, pp. e380
  * [26] John T Ormerod and Matt P Wand  “Gaussian variational approximate inference for generalized linear mixed models”  In _Journal of Computational and Graphical Statistics_ 21, 2012, pp. 2–17
  * [27] Danny Pfeffermann  “New important developments in small area estimation”  In _Statistical Science_ 28, 2013, pp. 40–68
  * [28] NG Narasimha Prasad and Jon NK Rao  “The estimation of the mean squared error of small-area estimators”  In _Journal of the American statistical association_ 85 Taylor & Francis, 1990, pp. 163–171
  * [29] Sara Geer and Patric Müller  “Quasi-likelihood and/or robust estimation in high dimensions”  In _Statistical Science_ 27, 2012, pp. 469–480
  * [30] Edward F Vonesh, Hao Wang, Lei Nie and Dibyen Majumdar  “Conditional second-order generalized estimating equations for generalized linear and nonlinear mixed-effects models”  In _Journal of the American Statistical Association_ 97, 2002, pp. 271–283

## References

  * [31] Ishwar V Basawa and David John Scott  “Asymptotic optimal inference for non-ergodic models”  Springer Science & Business Media, 2012
  * [32] Douglas Bates, Martin Mächler, Ben Bolker and Steve Walker  “Fitting Linear Mixed-Effects Models Using lme4”  In _Journal of Statistical Software_ 67.1, 2015, pp. 1–48  DOI: 10.18637/jss.v067.i01
  * [33] JR Blum, H Chernoff, M Rosenblatt and H Teicher  “Central limit theorems for interchangeable processes”  In _Canadian Journal of Mathematics_ 10 Cambridge University Press, 1958, pp. 222–229
  * [34] Norman E Breslow and David G Clayton  “Approximate inference in generalized linear mixed models”  In _Journal of the American statistical Association_ 88, 1993, pp. 9–25
  * [35] Norman E Breslow and Xihong Lin  “Bias correction in generalised linear mixed models with a single component of dispersion”  In _Biometrika_ 82, 1995, pp. 81–91
  * [36] Mollie E. Brooks et al.  “glmmTMB Balances Speed and Flexibility Among Packages for Zero-inflated Generalized Linear Mixed Modeling”  In _The R Journal_ 9, 2017, pp. 378–400
  * [37] Lucien Le Cam and Grace L. Yang  “On the Preservation of Local Asymptotic Normality under Information Loss”  In _The Annals of Statistics_ 16 Institute of Mathematical Statistics, 1988, pp. 483–520
  * [38] Jing Cheng et al.  “Real longitudinal data analysis for real people: building a good enough mixed model”  In _Statistics in medicine_ 29, 2010, pp. 504–520
  * [39] Peter J Downey  “Distribution-free bounds on the expectation of the maximum with scheduling applications”  In _Operations Research Letters_ 9, 1990, pp. 189–201
  * [40] Yingying Fan and Runze Li  “Variable selection in linear mixed effects models”  In _Annals of statistics_ 40, 2012, pp. 2043–2068
  * [41] Charles R Henderson  “Sire evaluation and genetic trends”  In _Journal of Animal Science_ 1973, 1973, pp. 10–41
  * [42] Francis K C Hui  “On the use of a penalized quasilikelihood information criterion for generalized linear mixed models”  In _Biometrika_ 108, 2020, pp. 353–365
  * [43] Francis K C Hui, Samuel Müller and A H Welsh  “Joint selection in mixed models using regularized PQL”  In _Journal of the American Statistical Association_ 112, 2017, pp. 1323–1333
  * [44] Francis K C Hui, Samuel Müller and Alan H Welsh  “Random effects misspecification can have severe consequences for random effects inference in linear mixed models”  In _International Statistical Review_ 89, 2021, pp. 186–206
  * [45] Jiming Jiang  “Empirical best prediction for small-area inference based on generalized linear mixed models”  In _Journal of Statistical Planning and Inference_ 111.1-2 Elsevier, 2003, pp. 117–127
  * [46] Jiming Jiang, Haomiao Jia and Hegang Chen  “MAXIMUM POSTERIOR ESTIMATION OF RANDOM EFFECTS IN GENERALIZED LINEAR MIXED MODELS”  In _Statistica Sinica_ 11.1 Institute of Statistical Science, Academia Sinica, 2001, pp. 97–120
  * [47] Jiming Jiang, Matt P. Wand and Aishwarya Bhaskaran  “Usable and precise asymptotics for generalized linear mixed model analysis and design”  In _Journal of the Royal Statistical Society: Series B_ 84, 2022, pp. 55–82
  * [48] Raghu N Kackar and David A Harville  “Approximations for standard errors of estimators of fixed and random effects in mixed linear models”  In _Journal of the American Statistical Association_ 79 Taylor & Francis, 1984, pp. 853–862
  * [49] Łukasz Kidziński, Francis KC Hui, David I Warton and Trevor J Hastie  “Generalized Matrix Factorization: efficient algorithms for fitting generalized linear latent variable models to large data arrays”  In _The Journal of Machine Learning Research_ 23, 2022, pp. 13211–13239
  * [50] Ziyang Lyu and A H Welsh  “Asymptotics for EBLUPs: Nested Error Regression Models”  In _Journal of the American Statistical Association_ 117, 2021, pp. 1–15
  * [51] Ziyang Lyu and A H Welsh  “Increasing Cluster Size Asymptotics for Nested Error Regression Models”  In _Journal of Statistical Planning and Inference_ 217, 2021, pp. 52–68
  * [52] Charles E McCulloch and Shayle R Searle  “Generalized, linear, and mixed models”  John Wiley & Sons, 2004
  * [53] Lei Nie  “Convergence rate of MLE in generalized linear and nonlinear mixed-effects models: Theory and applications”  In _Journal of Statistical Planning and Inference_ 137, 2007, pp. 1787–1804
  * [54] HE Ogden  “On asymptotic validity of naive inference with an approximate likelihood”  In _Biometrika_ 104, 2017, pp. 153–164
  * [55] Helen Ogden  “On the error in Laplace approximations of high-dimensional integrals”  In _Stat_ 10, 2021, pp. e380
  * [56] John T Ormerod and Matt P Wand  “Gaussian variational approximate inference for generalized linear mixed models”  In _Journal of Computational and Graphical Statistics_ 21, 2012, pp. 2–17
  * [57] Danny Pfeffermann  “New important developments in small area estimation”  In _Statistical Science_ 28, 2013, pp. 40–68
  * [58] NG Narasimha Prasad and Jon NK Rao  “The estimation of the mean squared error of small-area estimators”  In _Journal of the American statistical association_ 85 Taylor & Francis, 1990, pp. 163–171
  * [59] Sara Geer and Patric Müller  “Quasi-likelihood and/or robust estimation in high dimensions”  In _Statistical Science_ 27, 2012, pp. 469–480
  * [60] Edward F Vonesh, Hao Wang, Lei Nie and Dibyen Majumdar  “Conditional second-order generalized estimating equations for generalized linear and nonlinear mixed-effects models”  In _Journal of the American Statistical Association_ 97, 2002, pp. 271–283

Supplementary Material

In the developments, we prove all results below assuming the working
dispersion parameter $\hat{\phi}$ is equal to the true dispersion parameter
$\dot{\phi}$. Then for the general result using any $O_{p}(1)$ working
$\hat{\phi}$, we note that solving

$\displaystyle\nabla Q(\hat{\bm{\theta}})$
$\displaystyle=\left[\begin{matrix}{\hat{\phi}}^{-1}\bm{X}^{\top}\\{\bm{y}-\bm{\mu}(\hat{\bm{\theta}})\\}\\\
{\hat{\phi}}^{-1}\bm{Z}^{\top}\\{\bm{y}-\bm{\mu}(\hat{\bm{\theta}})\\}-(\bm{I}_{m}\otimes\hat{\bm{G}}^{-1})\hat{\bm{b}}\end{matrix}\right]=\bm{0}_{(m+1)p}$

for $\hat{\bm{\theta}}$ is equivalent to solving

$\displaystyle\nabla Q(\hat{\bm{\theta}})$
$\displaystyle=\left[\begin{matrix}\dot{\phi}^{-1}\bm{X}^{\top}\\{\bm{y}-\bm{\mu}(\hat{\bm{\theta}})\\}\\\
\dot{\phi}^{-1}\bm{Z}^{\top}\\{\bm{y}-\bm{\mu}(\hat{\bm{\theta}})\\}-(\bm{I}_{m}\otimes\hat{\bm{G}}_{s}^{-1})\hat{\bm{b}}\end{matrix}\right]=\bm{0}_{(m+1)p},$

where $\hat{\bm{G}}_{s}=\dot{\phi}\hat{\phi}^{-1}\hat{\bm{G}}$, whose inverse
is still $O_{p}(1)$ and positive definite. This is equivalent to setting
$\hat{\phi}$ to $\dot{\phi}$ and scaling $\hat{\bm{G}}$ by
$\dot{\phi}\hat{\phi}^{-1}$. The general result then follows since the results
proved under $\hat{\phi}=\dot{\phi}$ hold for any $\hat{\bm{G}}$ that has an
$O_{p}(1)$, positive definite inverse.

###  S0.1 Bias and Identifiability in the Conditional Regime

By differentiating (2), we see that the PQL estimators satisfy
$\sum_{i=1}^{m}{\hat{\phi}}^{-1}\bm{X}_{i}^{\top}\\{\bm{y}_{i}-\bm{\mu}_{i}(\hat{\bm{\theta}})\\}=\bm{0}$
and
${\hat{\phi}}^{-1}\bm{Z}_{i}^{\top}\\{\bm{y}_{i}-\bm{\mu}_{i}(\hat{\bm{\theta}})\\}-\bm{G}^{-1}\hat{\bm{b}}_{i}=\bm{0}$,
$i=1,\ldots,m$. Summing both sides of the second equation across all $i$,
since $\bm{X}_{i}=\bm{Z}_{i}$, it follows that
$\sum_{i=1}^{m}\hat{\bm{b}}_{i}=\bm{0}_{p}$. That is, the PQL estimators of
the random effects must satisfy a sum-to-zero constraint regardless of the
underlying true parameter values. Under a conditional regime, this induces an
asymptotic bias as captured by the term
$\bm{1}^{*}_{m}\otimes(m^{-1}\sum_{i=1}^{m}\dot{\bm{b}}_{i})$ in Theorem 1,
which can be interpreted as shifting the mean of the random effects into the
corresponding fixed effects. We can deal with the bias by reparametrising the
model _a priori_ to satisfy a sum-to-zero constraint. That is, we can define a
reparametrized vector of true values $\dot{\bm{\theta}}^{*}$ which satisfy
$\bm{1}^{*}_{m}\otimes(m^{-1}\sum_{i=1}^{m}\dot{\bm{b}}_{i}^{*})=\bm{0}_{(m+1)p}$,
and the PQL estimator will then be asymptotically normally distributed
centered around $\dot{\bm{\theta}}^{*}$. Furthermore, Theorem 1 remains
practically useful as, for any given sample size, we can always reparameterise
the GLMM to satisfy this identifiability constraint.

The asymptotic bias discussed above is analogous to that seen in a over-
parametrized one-way analysis of variance (ANOVA) model. That is, in the ANOVA
model one can always reparametrise to satisfy a sum-to-zero constraint, and
the corresponding estimator is consistent for this vector of the
reparametrized true values. Note however that when we work unconditionally
(Section 4), reparametrising in this way will lead to a different model to the
original, since the clusters are no longer independent.

##  S1 Proofs for Consistency

To establish our large sample distributional results, we first require the
following consistency result.

######  Lemma 1.

Suppose Conditions (C1)-(C5) hold and $mn_{L}^{-2}\rightarrow 0$. Then, as
$m,n_{L}\rightarrow\infty$ and unconditional on the random effects
$\dot{\bm{b}}$, $\|\hat{\bm{\theta}}-\dot{\bm{\theta}}\|_{\infty}=o_{p}(1)$.

These results are required to control the remainder term in the Taylor
expansions we use to derive the distributional results in Section S2. To prove
the result, we wish to show that for any given $\epsilon>0$, there exists a
large enough constant $C>0$ such that, for large $m,n_{L}$, we have

$\displaystyle
P\left\\{\underset{\|\bm{u}\|_{\infty}=C}{\sup}Q(\dot{\bm{\theta}}+\delta_{m,n_{L}}^{-1}\bm{u})<Q(\dot{\bm{\theta}})\right\\}\geq
1-\epsilon,$

for some positive, unbounded, monotonically increasing sequence
$\delta_{m,n_{L}}$. The above result implies that with probability tending to
one, there exists a local maximum $\hat{\bm{\theta}}$ in the ball
$\\{\dot{\bm{\theta}}+\delta_{m,n_{L}}^{-1}\bm{u}:\|\bm{u}\|_{\infty}\leq
C\\}$ so that
$\|\delta_{m,n_{L}}(\hat{\bm{\theta}}-\dot{\bm{\theta}})\|_{\infty}=O_{p}(1)$,
and thus $\|\hat{\bm{\theta}}-\dot{\bm{\theta}}\|_{\infty}=o_{p}(1)$.

Consider the difference $Q(\dot{\bm{\theta}}+\bm{u})-Q(\dot{\bm{\theta}})$. By
a Taylor expansion, we obtain

$\displaystyle Q(\dot{\bm{\theta}}+\bm{u})-Q(\dot{\bm{\theta}})$
$\displaystyle=\bm{u}^{\top}\\{\nabla
Q(\dot{\bm{\theta}})\\}-0.5\bm{u}^{\top}\\{-\nabla^{2}Q(\bar{\bm{\theta}})\\}\bm{u}.$
(S1.1)

where $\bar{\bm{\theta}}$ lies on the line segment joining $\dot{\bm{\theta}}$
and $\dot{\bm{\theta}}+\bm{u}$. If we can prove that (S1.1) is negative as
$m,n_{L}\rightarrow\infty$ for any choice of $C$, then there must exist some
$\delta_{m,n_{L}}$ such that
$Q(\dot{\bm{\theta}}+\delta_{m,n_{L}}^{-1}\bm{u})-Q(\dot{\bm{\theta}})$ is
negative for large enough $C$, and the required result follows. We have

$\displaystyle\nabla Q(\dot{\bm{\theta}})$
$\displaystyle=\left[\begin{matrix}\dot{\phi}^{-1}\bm{X}^{\top}(\bm{y}-\dot{\bm{\mu}})\\\
\dot{\phi}^{-1}\bm{Z}^{\top}(\bm{y}-\dot{\bm{\mu}})-(\bm{I}_{m}\otimes\hat{\bm{G}}^{-1})\dot{\bm{b}}\end{matrix}\right]=\left[\begin{matrix}\dot{\phi}^{-1}\bm{X}^{\top}(\bm{y}-\dot{\bm{\mu}})\\\
\dot{\phi}^{-1}\bm{Z}^{\top}(\bm{y}-\dot{\bm{\mu}})\end{matrix}\right]+\left[\begin{matrix}\bm{0}_{p}\\\
-(\bm{I}_{m}\otimes\hat{\bm{G}}^{-1})\dot{\bm{b}}\end{matrix}\right]$
$\displaystyle\triangleq\bm{\lambda}_{1}+\bm{\lambda}_{2},$

and

$\displaystyle-\nabla^{2}Q(\bar{\bm{\theta}})$
$\displaystyle=\left[\begin{matrix}\bm{X}^{\top}\bar{\bm{W}}\bm{X}&\bm{X}_{1}^{\top}\bar{\bm{W}}_{1}\bm{X}_{1}&\cdots&\bm{X}_{m}^{\top}\bar{\bm{W}}_{m}\bm{X}_{m}\\\
\bm{X}_{1}^{\top}\bar{\bm{W}}_{1}\bm{X}_{1}&\bm{X}_{1}^{\top}\bar{\bm{W}}_{1}\bm{X}_{1}+\hat{\bm{G}}^{-1}&&\bm{0}\\\
\vdots&\vdots&\ddots&\vdots\\\
\bm{X}_{m}^{\top}\bar{\bm{W}}_{m}\bm{X}_{m}&\bm{0}&&\bm{X}_{m}^{\top}\bar{\bm{W}}_{m}\bm{X}_{m}+\hat{\bm{G}}^{-1}\end{matrix}\right]$
$\displaystyle=\left[\begin{matrix}\bm{X}^{\top}\bar{\bm{W}}\bm{X}&\bm{X}_{1}^{\top}\bar{\bm{W}}_{1}\bm{X}_{1}&\cdots&\bm{X}_{m}^{\top}\bar{\bm{W}}_{m}\bm{X}_{m}\\\
\bm{X}_{1}^{\top}\bar{\bm{W}}_{1}\bm{X}_{1}&\bm{X}_{1}^{\top}\bar{\bm{W}}_{1}\bm{X}_{1}&&\bm{0}\\\
\vdots&\vdots&\ddots&\vdots\\\
\bm{X}_{m}^{\top}\bar{\bm{W}}_{m}\bm{X}_{m}&\bm{0}&&\bm{X}_{m}^{\top}\bar{\bm{W}}_{m}\bm{X}_{m}\end{matrix}\right]+\text{blockdiag}(\bm{0}_{p},\bm{I}_{m}\otimes\hat{\bm{G}}^{-1})$
$\displaystyle\triangleq\bm{\Gamma}_{1}(\bar{\bm{\theta}})+\bm{\Gamma}_{2},$

where
$\bar{\bm{W}}_{i}=\dot{\phi}^{-1}\mbox{diag}{\\{a^{\prime\prime}(\bar{\eta}_{i1}),\ldots,a^{\prime\prime}(\bar{\eta}_{in_{i}})\\}}$
and
$\bar{\bm{W}}=\dot{\phi}^{-1}\mbox{diag}{\\{a^{\prime\prime}(\bar{\eta}_{11}),\ldots,a^{\prime\prime}(\bar{\eta}_{mn_{m}})\\}}$.
Also, let $\bm{\Gamma}_{1}(\dot{\bm{\theta}})+\bm{\Gamma}_{2}$ denote the
analogous decomposition of $-\nabla^{2}Q(\dot{\bm{\theta}})$. For both the
conditional and unconditional regimes, we will prove that the second term is
positive and dominates the first. However, the treatment of the terms differs
between the two cases, and as such the proofs will need to be dealt with
separately. In the following three sections, we will first treat the Poisson
pure random intercept example, followed by the more general conditional and
unconditional regimes.

Before proceeding, we demonstrate an inequality that is used in the proofs
below. Write $\bm{u}=(\bm{u}_{1}^{\top},\bm{u}_{2}^{\top})^{\top}$,
$\bm{u}_{2}=(\bm{u}_{21}^{\top},\ldots,\bm{u}_{2m}^{\top})^{\top}$. First, for
any $\bm{\theta}$ we have

$\displaystyle\bm{u}^{\top}\bm{\Gamma}_{1}(\bm{\theta})\bm{u}=\bm{u}_{1}^{\top}\bm{X}^{\top}\bm{W}\bm{X}\bm{u}_{1}+2\bm{u}_{1}^{\top}\bm{X}^{\top}\bm{W}\bm{Z}\bm{u}_{2}+\bm{u}_{2}^{\top}\bm{Z}^{\top}\bm{W}\bm{Z}\bm{u}_{2}^{\top}\geq
0.$

Next, we have

$\displaystyle\bm{u}^{\top}\bm{\Gamma}_{1}(\bar{\bm{\theta}})\bm{u}-c_{0}^{2}\bm{u}^{\top}\bm{\Gamma}_{1}(\bm{\theta})\bm{u}$
$\displaystyle=\bm{u}_{1}^{\top}\bm{X}^{\top}(\bar{\bm{W}}-c_{0}^{2}\bm{W})\bm{X}\bm{u}_{1}+2\bm{u}_{1}^{\top}\bm{X}^{\top}(\bar{\bm{W}}-c_{0}^{2}\bm{W})\bm{Z}\bm{u}_{2}+\bm{u}_{2}^{\top}\bm{Z}^{\top}(\bar{\bm{W}}-c_{0}^{2}\bm{W})\bm{Z}\bm{u}_{2}^{\top}.$

If we denote $\bm{W}^{*}=\bar{\bm{W}}-c_{0}^{2}\bm{W}$, then by Condition (C1)
$\bm{W}^{*}$ is a diagonal matrix with non-negative entries as the entries of
$c_{0}^{2}\bm{W}$ are upper bounded by the smallest component in
$\bar{\bm{W}}$. Therefore

$\displaystyle\bm{u}^{\top}\bm{\Gamma}_{1}(\bar{\bm{\theta}})\bm{u}-c_{0}^{2}\bm{u}^{\top}\bm{\Gamma}_{1}(\bm{\theta})\bm{u}=\bm{u}_{1}^{\top}\bm{X}^{\top}\bm{W}^{*}\bm{X}\bm{u}_{1}+2\bm{u}_{1}^{\top}\bm{X}^{\top}\bm{W}^{*}\bm{Z}\bm{u}_{2}+\bm{u}_{2}^{\top}\bm{Z}^{\top}\bm{W}^{*}\bm{Z}\bm{u}_{2}^{\top}\geq
0,$

so that $\bm{u}^{\top}\bm{\Gamma}_{1}(\bar{\bm{\theta}})\bm{u}\geq
c_{0}^{2}\bm{u}^{\top}\bm{\Gamma}_{1}(\bm{\theta})\bm{u}$. Finally, note that
we can choose $\bm{\theta}=\dot{\bm{\theta}}$ or
$\bm{\theta}=E(\dot{\bm{\theta}})$ without altering the above argument.

###  S1.1 Poisson pure random intercept example

We begin with the Poisson pure random intercept example, which gives insight
and covers a case where $\bm{X}_{i}\neq\bm{Z}_{i}$. The following result is
unconditional on the random effects $\dot{\bm{b}}$.

######  Lemma 2.

Assume Conditions (C1)-(C5) hold, and let $mn^{-2}\rightarrow 0$. Then for the
Poisson pure random intercept model, as $m,n\rightarrow\infty$ and
unconditional on the random effects $\dot{\bm{b}}$, it holds that
$\|\hat{\bm{\theta}}-\dot{\bm{\theta}}\|_{\infty}=o_{p}(1)$.

###### Proof.

Let $\bm{u}=\bm{u}_{2}=({u}_{21},\ldots,{u}_{2m})^{\top}$,
$\bm{\theta}=\bm{b}=(b_{1},\ldots,b_{m})^{\top}$,
$\hat{\bm{G}}=\hat{\sigma}_{b}^{2}$ (a scalar),
$-\nabla^{2}Q(\bar{\bm{\theta}})=\mbox{diag}{(ne^{\bar{b}_{1}}+\hat{\sigma}_{b}^{-2},\ldots,ne^{\bar{b}_{m}}+\hat{\sigma}_{b}^{-2})}\equiv\bm{\Gamma}_{1}(\bar{\bm{\theta}})+\bm{\Gamma}_{2}$,
and

$\displaystyle\nabla
Q(\dot{\bm{\theta}})=\left[\begin{matrix}\sum_{j=1}^{n}(y_{1j}-e^{\dot{b}_{1}})-\hat{\sigma}_{b}^{-2}\dot{b}_{1}\\\
\vdots\\\
\sum_{j=1}^{n}(y_{mj}-e^{\dot{b}_{m}})-\hat{\sigma}_{b}^{-2}\dot{b}_{m}\end{matrix}\right]\equiv\bm{\lambda}_{1}+\bm{\lambda}_{2}.$

Let $\bm{M}=E\\{\mbox{diag}{(ne^{\dot{b}_{1}},\ldots,ne^{\dot{b}_{m}})}\\}$.
Then
$\bm{M}=\text{Var}\\{\dot{\phi}^{-1}\bm{Z}^{\top}(\bm{y}-\dot{\bm{\mu}})\\}$.
By Condition (C1),
$c_{0}^{2}\bm{u}^{\top}\bm{M}\bm{u}\leq\bm{u}^{\top}\bm{\Gamma}_{1}(\bar{\bm{\theta}})\bm{u}$.
Next, let $\lambda=\dot{\hat{\sigma}}_{b}^{2}\hat{\sigma}_{b}^{-2}$. Then
$\text{Var}(\bm{\lambda}_{2})=\dot{\hat{\sigma}}_{b}^{2}\hat{\sigma}_{b}^{-4}\bm{I}_{m}$
and

$\displaystyle\lambda^{-1}\bm{u}_{2}^{\top}(\dot{\hat{\sigma}}_{b}^{2}\hat{\sigma}_{b}^{-4}\bm{I}_{m})\bm{u}_{2}=\bm{u}_{2}^{\top}(\hat{\sigma}_{b}^{-2}\bm{I}_{m})\bm{u}_{2}=\bm{u}^{\top}\bm{\Gamma}_{2}\bm{u}.$

Finally, by the laws of iterated expectation and variance we have

$\displaystyle E\\{\nabla Q(\dot{\bm{\theta}})\nabla
Q(\dot{\bm{\theta}})^{\top}\\}$ $\displaystyle=\text{Var}\\{\nabla
Q(\dot{\bm{\theta}})\\}$ $\displaystyle=E[\text{Var}\\{\nabla
Q(\dot{\bm{\theta}})|\dot{\bm{b}}\\}]+\text{Var}[E\\{\nabla
Q(\dot{\bm{\theta}})|\dot{\bm{b}}\\}]$
$\displaystyle=E\\{\text{Var}(\bm{\lambda}_{1}|\dot{\bm{b}})\\}+\text{Var}(\bm{\lambda}_{2})$
$\displaystyle=\text{Var}(\bm{\lambda}_{1})+\text{Var}(\bm{\lambda}_{2}).$

Therefore, we have that

$\displaystyle\bm{u}^{\top}\\{-\nabla^{2}Q(\bar{\bm{\theta}})\\}\bm{u}$
$\displaystyle\geq\min(\lambda^{-1},c_{0}^{2})\\{\bm{u}^{\top}\bm{M}\bm{u}+\bm{u}_{2}^{\top}(\dot{\hat{\sigma}}_{b}^{2}\hat{\sigma}_{b}^{-4}\bm{I}_{m})\bm{u}_{2}\\}$
$\displaystyle=\min(\lambda^{-1},c_{0}^{2})E\\{\nabla
Q(\dot{\bm{\theta}})\nabla Q(\dot{\bm{\theta}})^{\top}\\},$

where the latter and hence former term grows at the same rate as
$\\{\bm{u}^{\top}\nabla Q(\dot{\bm{\theta}})\\}^{2}$. Since at least one
component of $\bm{u}$ equals $\pm C$, for any given $\bm{u}$,
$\bm{u}^{\top}\\{-\nabla^{2}Q(\bar{\bm{\theta}})\\}\bm{u}$ is at least of
order $O_{p}(m)$ in probability and hence always dominates.

Since the choice of which $|u_{2i}|=C$ is arbitrary however, we also need to
make sure that the $m$th order statistic
$\underset{i\in\\{1,\ldots,m\\}}{\max}[\\{\sum_{j=1}^{n}(y_{ij}-e^{\dot{b}_{i}})-\hat{\sigma}_{b}^{-2}\dot{b}_{i}\\}/(ne^{\dot{b}_{i}}+\hat{\sigma}_{b}^{-2})]$,
which grows with the dimension, is of order $o_{p}(1)$. We know that the
leading term in (3b) is
$(\bm{Z}^{\top}\dot{\bm{W}}\bm{Z})^{-1}\\{\dot{\phi}^{-1}\bm{Z}^{\top}(\bm{y}-\dot{\bm{\mu}})\\}$
when $mn^{-1}\rightarrow\infty$; for this Poisson random intercept example, up
to some smaller order terms, this simplifies to the ratio
$\\{\sum_{j=1}^{n}(y_{ij}-e^{\dot{b}_{i}})-\hat{\sigma}_{b}^{-2}\dot{b}_{i}\\}/(ne^{\dot{b}_{i}}+\hat{\sigma}_{b}^{-2})$.
Intuitively then, proving a result for
$\|\hat{\bm{\theta}}-\dot{\bm{\theta}}\|_{\infty}$ should involve studying
$\underset{i\in\\{1,\ldots,m\\}}{\max}[\\{\sum_{j=1}^{n}(y_{ij}-e^{\dot{b}_{i}})-\hat{\sigma}_{b}^{-2}\dot{b}_{i}\\}/(ne^{\dot{b}_{i}}+\hat{\sigma}_{b}^{-2})]$.

Put another way, consider the set of $\bm{u}$ such that one component of
$\bm{u}$ equals $\pm C$ and zero elsewhere. When $C$ is the $i$th component of
$\bm{u}$, this corresponds to deviating away from $\dot{\bm{\theta}}$ in the
$i$th direction. In this case, we need
$C\\{\sum_{j=1}^{n}(y_{ij}-e^{\dot{b}_{i}})-\hat{\sigma}_{b}^{-2}\dot{b}_{i}\\}$
to be dominated by $C^{2}ne^{\dot{b}_{i}}$ for any $C$ and all $m,n$ large
enough, i.e.,
$\\{\sum_{j=1}^{n}(y_{ij}-e^{\dot{b}_{i}})-\hat{\sigma}_{b}^{-2}\dot{b}_{i}\\}/ne^{\dot{b}_{i}}=o_{p}(1)$.
This is indeed true as this ratio is $O_{p}(n^{-1/2})$, since
$\sum_{j=1}^{n}(y_{ij}-e^{\dot{b}_{i}})-\hat{\sigma}_{b}^{-2}\dot{b}_{i}=O_{p}(n^{1/2})$
due to conditional independence and Chebyshev’s inequality, and
$e^{\dot{b}_{i}}=O_{p}(1)$. However, although the ratio is of order
$O_{p}(n^{-1/2})$, for any given $m,n$ there is still a positive probability
that the ratio (a random variable) is greater than one in magnitude. On the
other hand, for the consistency argument to hold we need to make sure the
ratio is smaller than one in magnitude for all $m$ directions with probability
tending to one, as $m,n\rightarrow\infty$. In particular, it is sufficient for
the maximum of $m$ of these ratios to be $o_{p}(1)$: this maximum grows with
$m$, corresponding to the number of directions we need to bound. Intuitively,
this should hold if $m$ does not grow too fast relative to $n$.

Now, [39] proves that the maximum over $m$ realisations of independently and
identically distributed random variables with a finite $q$th moment is
$o_{p}(m^{1/q})$. By Condition (C5), the ratio
$n^{1/2}\\{\sum_{j=1}^{n}(y_{ij}-e^{\dot{b}_{i}})-\hat{\sigma}_{b}^{-2}\dot{b}_{i}\\}/(ne^{\dot{b}_{i}}+\hat{\sigma}_{b}^{-2})$
has finite fourth moments for all $i$ and $n$. Thus, the maximum of these
(normalised) ratios over $m$ clusters is of order $o_{p}(m^{1/4})$. As a
result, the maximum ratio of interest is $o_{p}(m^{1/4}n^{-1/2})$. Therefore,
when $mn^{-2}\rightarrow 0$, there exists $\delta_{m,n}$ such that we can
always choose a large enough $C$ for $\delta_{m,n}^{-1}\bm{u}^{\top}\nabla
Q(\dot{\bm{\theta}})$ to be dominated by
$\delta_{m,n}^{-2}\bm{u}^{\top}\\{-\nabla^{2}Q(\bar{\bm{\theta}})\\}\bm{u}$,
and hence
$\|\delta_{m,n}(\hat{\bm{\theta}}-\dot{\bm{\theta}})\|_{\infty}=O_{p}(1)$ as
required.

∎

To conclude this section, we remark that although $mn^{-2}\rightarrow 0$ is
needed for the consistency and thus distributional result, this is a
sufficient condition. Intuitively, in the Poisson pure random effects model
there are no fixed parameters to estimate, and the estimate of the random
effects for each cluster only depends on observations in that cluster. Thus,
the relative rates of $m$ and $n$ should not matter for a distributional
result concerning a finite subset of the random effects.

###  S1.2 Conditional on the Random Effects

In this section, we prove the consistency result under the conditional regime.
In the conditional regime, we assume without loss of generality that
$\sum_{i=1}^{m}\dot{\bm{b}}_{i}=\bm{0}_{p}$, recalling that we can always
reparametrise the random effect coefficients so this holds.

Let $\bm{M}=\bm{\Gamma}_{1}(\dot{\bm{\theta}})$. Then
$\bm{M}=\text{Var}(\bm{\lambda}_{1}|\dot{\bm{b}})=E(\bm{\lambda}_{1}\bm{\lambda}_{1}^{\top}|\dot{\bm{b}})$
since $E(\bm{\lambda}_{1}|\dot{\bm{b}})=\bm{0}_{(m+1)p}$. By Condition (C1),
we have
$c_{0}^{2}\bm{u}^{\top}\bm{M}\bm{u}\leq\bm{u}^{\top}\bm{\Gamma}_{1}(\bar{\bm{\theta}})\bm{u}$.

We now consider two cases: the special case when $\bm{u}_{1}=-\bm{u}_{2i}$ for
all $i$, and when this is not the case. For the former, we have
$\bm{u}^{\top}\bm{\lambda}_{1}=\bm{u}^{\top}\bm{M}\bm{u}=0$. Then we must
examine $\bm{u}^{\top}\bm{\lambda}_{2}$ and
$\bm{u}^{\top}\bm{\Gamma}_{2}\bm{u}$. In this case, we have
$\bm{u}^{\top}\bm{\lambda}_{2}=\sum_{i=1}^{m}\bm{u}_{2i}^{\top}\hat{\bm{G}}^{-1}\dot{\bm{b}}_{i}=-\bm{u}_{1}^{\top}\hat{\bm{G}}^{-1}\sum_{i=1}^{m}\dot{\bm{b}}_{i}=0$,
and
$\bm{u}^{\top}\bm{\Gamma}_{2}\bm{u}=m\bm{u}_{1}^{\top}\hat{\bm{G}}^{-1}\bm{u}_{1}>0$
since $\hat{\bm{G}}$ is a positive definite matrix. Thus the difference (S1.1)
is negative for large enough $m,n_{L}$ and any choice of constant $C$.

Next, consider the case when $\bm{u}_{1}=-\bm{u}_{2i}$ for all $i$ does not
hold. Under this setting, as $\bm{\Gamma}_{2}$ is a positive semi-definite
matrix, we still have
$\bm{u}^{\top}\\{-\nabla^{2}Q(\bar{\bm{\theta}})\\}\bm{u}\geq\bm{u}^{\top}\bm{\Gamma}_{1}(\bar{\bm{\theta}})\bm{u}\geq
c_{0}^{2}\bm{u}^{\top}\bm{M}\bm{u}$, where the last and hence former terms
grow at the same rate as $(\bm{u}^{\top}\bm{\lambda}_{1})^{2}$. Since at least
one component of $\bm{u}$ equals $\pm C$, by Conditions (C1)-(C3) we have that
$\bm{u}^{\top}\\{-\nabla^{2}Q(\bar{\bm{\theta}})\\}\bm{u}$ is at least of
order $O_{p}(n_{L})$, and always dominates since
$\bm{u}^{\top}\bm{\lambda}_{2}=O_{p}(m)$ at most.

Since the choice of $\bm{u}$ is arbitrary, we must take into account the
growth rate of the $m$th order statistic. That is, for any $1\leq k\leq p$, we
need
$\underset{i\in\\{1,\ldots,m\\}}{\max}[(\bm{X}_{i}^{\top}\dot{\bm{W}}_{i}\bm{X}_{i}+\hat{\bm{G}}^{-1})^{-1}\\{\dot{\phi}^{-1}\bm{X}_{i}^{\top}(\bm{y}_{i}-\dot{\bm{\mu}}_{i})-\hat{\bm{G}}^{-1}\dot{\bm{b}}_{i}\\}]_{[k]}=o_{p}(1)$,
as per the argument for the Poisson pure random intercept model. Since the
responses $y_{ij}$ are from the exponential family and thus the moment
generating function always exists, the maximum is of order
$o_{p}(m^{1/r}n_{L}^{-1/2})$ for any $r\in\mathbb{N}$ [39], and hence
$o_{p}(1)$ since $mn_{L}^{-1}\rightarrow 0$ by taking $r=2$, for example. Note
that the first $p$ components of $\nabla Q(\dot{\bm{\theta}})$, which are
associated with the fixed effects, do not need to be bounded in this way
because the dimension is fixed.

###  S1.3 Unconditional on the Random Effects

In this section, we prove the consistency result under the unconditional
regime. The main differences to the derivation under the conditional regime
arise from the treatment of $\bm{\lambda}_{2}$, and the distribution of
$\bm{y}$. In the unconditional regime it holds that
$\sum_{i=1}^{m}\dot{\bm{b}}_{i}=O_{p}(m^{1/2})$, while in the conditional
regime we impose a sum to zero constraint. Furthermore, in the unconditional
regime we bound $\bm{u}^{\top}\bm{\lambda}_{2}$ using its variance, while in
the conditional regime this is not possible because $\bm{\lambda}_{2}$ is not
a random variable. Finally, in the unconditional regime we cannot use the
properties of the exponential family to bound the $m$th order statistic,
instead requiring Condition (C5).

Let $\bm{M}=E\\{\bm{\Gamma}_{1}(\dot{\bm{\theta}})\\}$. Then
$\bm{M}=\text{Var}(\bm{\lambda}_{1})=E(\bm{\lambda}_{1}\bm{\lambda}_{1}^{\top})$
since $E(\bm{\lambda}_{1})=\bm{0}_{(m+1)p}$. By Condition (C1),
$c_{0}^{2}\bm{u}^{\top}\bm{M}\bm{u}\leq\bm{u}^{\top}\bm{\Gamma}_{1}(\bar{\bm{\theta}})\bm{u}$.

We consider two cases: the special case when $\bm{u}_{1}=-\bm{u}_{2i}$ for all
$i$, and when this is not the case. In the former, we have
$\bm{u}^{\top}\bm{\lambda}_{1}=\bm{u}^{\top}\bm{M}\bm{u}=0$. Thus we must
examine $\bm{u}^{\top}\bm{\lambda}_{2}$ and
$\bm{u}^{\top}\bm{\Gamma}_{2}\bm{u}$. In this case, we have
$\bm{u}^{\top}\bm{\lambda}_{2}=\sum_{i=1}^{m}\bm{u}_{2i}^{\top}\hat{\bm{G}}^{-1}\dot{\bm{b}}_{i}=-\bm{u}_{1}^{\top}\hat{\bm{G}}^{-1}\sum_{i=1}^{m}\dot{\bm{b}}_{i}=O_{p}(m^{1/2})$,
and
$\bm{u}^{\top}\bm{\Gamma}_{2}\bm{u}=m\bm{u}_{1}^{\top}\hat{\bm{G}}^{-1}\bm{u}_{1}>0$
since $\hat{\bm{G}}$ is a positive definite matrix. Hence the difference
(S1.1) is negative for large enough $m,n_{L}$, and any choice of constant $C$.

Next, consider the case when $\bm{u}_{1}=-\bm{u}_{2i}$ for all $i$ does not
hold. Then we still have
$\bm{u}^{\top}\\{-\nabla^{2}Q(\bar{\bm{\theta}})\\}\bm{u}\geq
c_{0}^{2}\bm{u}^{\top}\bm{M}\bm{u}$. Letting
$\lambda=\lambda_{\max}(\hat{\bm{G}}^{-1}\dot{\bm{G}}\hat{\bm{G}}^{-1})/\lambda_{\min}(\hat{\bm{G}}^{-1})$,
we have

$\displaystyle\text{Var}(\bm{\lambda}_{2})=\bm{I}_{m}\otimes\hat{\bm{G}}^{-1}\dot{\bm{G}}\hat{\bm{G}}^{-1}$

and

$\displaystyle\lambda^{-1}\bm{u}_{2}^{\top}(\bm{I}_{m}\otimes\hat{\bm{G}}^{-1}\dot{\bm{G}}\hat{\bm{G}}^{-1})\bm{u}_{2}\leq\bm{u}_{2}^{\top}(\bm{I}_{m}\otimes\hat{\bm{G}}^{-1})\bm{u}_{2}=\bm{u}^{\top}\bm{\Gamma}_{2}\bm{u}.$

Now, by the laws of iterated expectation and variance,

$\displaystyle E\\{\nabla Q(\dot{\bm{\theta}})\nabla
Q(\dot{\bm{\theta}})^{\top}\\}$ $\displaystyle=\text{Var}\\{\nabla
Q(\dot{\bm{\theta}})\\}$ $\displaystyle=E[\text{Var}\\{\nabla
Q(\dot{\bm{\theta}})|\dot{\bm{b}}\\}]+\text{Var}[E\\{\nabla
Q(\dot{\bm{\theta}})|\dot{\bm{b}}\\}]$
$\displaystyle=E\\{\text{Var}(\bm{\lambda}_{1}|\dot{\bm{b}})\\}+\text{Var}(\bm{\lambda}_{2})$
$\displaystyle=\text{Var}(\bm{\lambda}_{1})+\text{Var}(\bm{\lambda}_{2}).$

Thus we have that

$\displaystyle\bm{u}^{\top}\\{-\nabla^{2}Q(\bar{\bm{\theta}})\\}\bm{u}$
$\displaystyle\geq\min(\lambda^{-1},c_{0}^{2})\\{\bm{u}^{\top}\bm{M}\bm{u}+\bm{u}_{2}^{\top}(\bm{I}_{m}\otimes\hat{\bm{G}}^{-1}\dot{\bm{G}}\hat{\bm{G}}^{-1})\bm{u}_{2}\\}$
$\displaystyle=\min(\lambda^{-1},c_{0}^{2})\bm{u}^{\top}E\\{\nabla
Q(\dot{\bm{\theta}})\nabla Q(\dot{\bm{\theta}})^{\top}\\}\bm{u},$

where the latter and hence former term grows at the same rate as
$\\{\bm{u}^{\top}\nabla Q(\dot{\bm{\theta}})\\}^{2}$. Since at least one
component of $\bm{u}$ equals $\pm C$, for any given $\bm{u}$ we have that
$\bm{u}^{\top}\\{-\nabla^{2}Q(\bar{\bm{\theta}})\\}\bm{u}$ is at least of
order $O_{p}(n_{L})$ and always dominates.

Since the choice of $\bm{u}$ is arbitrary, we must take into account the
growth rate of the $n$th order statistic. That is, for any $1\leq k\leq p$, we
require
$\underset{i\in\\{1,\ldots,m\\}}{\max}[(\bm{X}_{i}^{\top}\dot{\bm{W}}_{i}\bm{X}_{i}+\hat{\bm{G}}^{-1})^{-1}\\{\dot{\phi}^{-1}\bm{X}_{i}^{\top}(\bm{y}_{i}-\dot{\bm{\mu}}_{i})-\hat{\bm{G}}^{-1}\dot{\bm{b}}_{i}\\}]_{[k]}=o_{p}(1)$,
as per the argument for the Poisson pure random intercept model. By Condition
(C5), this term is of order $o_{p}(m^{1/4}n_{L}^{-1/2})$, and hence the result
follows. Note that the first $p$ components of $\nabla Q(\dot{\bm{\theta}})$,
which are associated with the fixed effects, do not need to be bounded in this
way because the dimension is fixed.

##  S2 Proofs of Distributional Results

For both the conditional and unconditional regimes, our proof relies on
examining the behaviour of the leading term in the Taylor expansion of the
estimating function. Under Conditions (C1) and (C3), we take the Taylor
expansion of $\nabla Q(\hat{\bm{\theta}})$ around $\dot{\bm{\theta}}$ and
obtain, as $m,n_{L}\rightarrow\infty$,

$\displaystyle\nabla Q(\hat{\bm{\theta}})=\bm{0}_{(m+1)p}=\nabla
Q(\dot{\bm{\theta}})+\nabla^{2}Q(\dot{\bm{\theta}})(\hat{\bm{\theta}}-\dot{\bm{\theta}})+\frac{1}{2}\bm{R}(\tilde{\bm{\theta}}),$
(S2.1)

where $\tilde{\bm{\theta}}$ is a $(m+1)p\times(m+1)p$ matrix with each row
lying on the line segment between $\dot{\bm{\theta}}$ and $\hat{\bm{\theta}}$
and $\bm{R}(\tilde{\bm{\theta}})$ is the remainder term. Rearranging, we have

$\displaystyle\hat{\bm{\theta}}-\dot{\bm{\theta}}$
$\displaystyle=-\\{\nabla^{2}Q(\dot{\bm{\theta}})\\}^{-1}\nabla
Q(\dot{\bm{\theta}})-\frac{1}{2}\\{\nabla^{2}Q(\dot{\bm{\theta}})\\}^{-1}\bm{R}(\tilde{\bm{\theta}}).$
(S2.2)

We show in Section S3 that the remainder term is of smaller order than the
leading term and thus negligible in the limit, in both the conditional and
unconditional regimes.

From (S2.2), to study the asymptotic behaviour of the PQL estimator we will
first apply the blockwise matrix inversion formula to obtain an expression for
$-\\{\nabla^{2}Q(\dot{\bm{\theta}})\\}^{-1}$. Using this result, we will then
obtain an expression for $-\\{\nabla^{2}Q(\dot{\bm{\theta}})\\}^{-1}\nabla
Q(\dot{\bm{\theta}})$, and subsequently study the asymptotic behaviour of each
constituent term. Note that since $\nabla Q(\dot{\bm{\theta}})$ is a
$(m+1)p$-vector and $-\\{\nabla^{2}Q(\dot{\bm{\theta}})\\}^{-1}$ is a
$(m+1)p\times(m+1)p$ matrix, we cannot simply take their limits as per
standard fixed dimension asymptotics. Instead, we must evaluate
$-\\{\nabla^{2}Q(\dot{\bm{\theta}})\\}^{-1}\nabla Q(\dot{\bm{\theta}})$ as a
whole.

We can write

$\displaystyle\nabla Q(\dot{\bm{\theta}})$
$\displaystyle=\left[\begin{matrix}\dot{\phi}^{-1}\bm{X}^{\top}(\bm{y}-\dot{\bm{\mu}})\\\
\dot{\phi}^{-1}\bm{Z}^{\top}(\bm{y}-\dot{\bm{\mu}})-(\bm{I}_{m}\otimes\hat{\bm{G}}^{-1})\dot{\bm{b}}\end{matrix}\right]=\left[\begin{matrix}\dot{\phi}^{-1}\sum_{i=1}^{m}\sum_{j=1}^{n_{i}}\bm{x}_{ij}(y_{ij}-\dot{\mu}_{ij})\\\
\dot{\phi}^{-1}\sum_{j=1}^{n_{1}}\bm{x}_{1j}(y_{1j}-\dot{\mu}_{1j})-\hat{\bm{G}}^{-1}\dot{\bm{b}}_{1}\\\
\vdots\\\
\dot{\phi}^{-1}\sum_{j=1}^{n_{m}}\bm{x}_{mj}(y_{mj}-\dot{\mu}_{mj})-\hat{\bm{G}}^{-1}\dot{\bm{b}}_{m}\end{matrix}\right]$
$\displaystyle\triangleq\left[\begin{matrix}\bm{S}_{1}\\\
\bm{S}_{21}+\bm{S}_{31}\\\ \vdots\\\
\bm{S}_{2m}+\bm{S}_{3m}\end{matrix}\right]\triangleq\left[\begin{matrix}\bm{S}_{1}\\\
\bm{S}_{4}+\bm{S}_{5}\end{matrix}\right]\triangleq\left[\begin{matrix}\bm{S}_{1}\\\
\bm{S}_{6}\end{matrix}\right],$ $\displaystyle\bm{B}(\dot{\bm{\theta}})$
$\displaystyle=-\nabla^{2}Q(\dot{\bm{\theta}})=\left[\begin{matrix}\bm{X}^{\top}\dot{\bm{W}}\bm{X}&\bm{X}^{\top}\dot{\bm{W}}\bm{Z}\\\
\bm{Z}^{\top}\dot{\bm{W}}\bm{X}&\bm{Z}^{\top}\dot{\bm{W}}\bm{Z}+\bm{I}_{m}\otimes\hat{\bm{G}}^{-1}\end{matrix}\right]\triangleq\left[\begin{matrix}\underset{p\times
p}{\bm{B}_{1}}&\underset{p\times mp}{\bm{B}_{2}}\\\ \underset{mp\times
p}{\bm{B}_{2}^{\top}}&\underset{mp\times
mp}{\bm{B}_{3}+\bm{B}_{4}}\end{matrix}\right].$

Letting
$\bm{C}=\bm{B}_{1}-\bm{B}_{2}(\bm{B}_{3}+\bm{B}_{4})^{-1}\bm{B}_{2}^{\top}$,
by the matrix block inversion formula we have

$\displaystyle\bm{B}^{-1}=\left[\begin{matrix}\bm{C}^{-1}&-\bm{C}^{-1}\bm{B}_{2}(\bm{B}_{3}+\bm{B}_{4})^{-1}\\\
-(\bm{B}_{3}+\bm{B}_{4})^{-1}\bm{B}_{2}^{\top}\bm{C}^{-1}&(\bm{B}_{3}+\bm{B}_{4})^{-1}+(\bm{B}_{3}+\bm{B}_{4})^{-1}\bm{B}_{2}^{\top}\bm{C}^{-1}\bm{B}_{2}(\bm{B}_{3}+\bm{B}_{4})^{-1}\end{matrix}\right].$
(S2.3)

Next, based on the forms of $\bm{B}_{2}$ and $(\bm{B}_{3}+\bm{B}_{4})$, we
obtain

$\displaystyle\bm{B}_{2}(\bm{B}_{3}+\bm{B}_{4})^{-1}=[\bm{I}_{p}-\hat{\bm{G}}^{-1}(\bm{X}^{\top}_{1}\dot{\bm{W}}_{1}\bm{X}_{1}+\hat{\bm{G}}^{-1})^{-1},\ldots,\bm{I}_{p}-\hat{\bm{G}}^{-1}(\bm{X}^{\top}_{m}\dot{\bm{W}}_{m}\bm{X}_{m}+\hat{\bm{G}}^{-1})^{-1}].$

Then since $\bm{Z}_{i}=\bm{X}_{i}$ for all $i$, we can show that

$\displaystyle\bm{B}_{2}(\bm{B}_{3}+\bm{B}_{4})^{-1}\bm{B}_{2}^{\top}$
$\displaystyle=\sum_{i=1}^{m}\bm{X}^{\top}_{i}\dot{\bm{W}}_{i}\bm{X}_{i}(\bm{X}^{\top}_{i}\dot{\bm{W}}_{i}\bm{X}_{i}+\hat{\bm{G}}^{-1})^{-1}\bm{X}^{\top}_{i}\dot{\bm{W}}_{i}\bm{X}_{i}$
$\displaystyle=\sum_{i=1}^{m}(\bm{X}^{\top}_{i}\dot{\bm{W}}_{i}\bm{X}_{i}+\hat{\bm{G}}^{-1}-\hat{\bm{G}}^{-1})(\bm{X}^{\top}_{i}\dot{\bm{W}}_{i}\bm{X}_{i}+\hat{\bm{G}}^{-1})^{-1}\bm{X}^{\top}_{i}\dot{\bm{W}}_{i}\bm{X}_{i}$
$\displaystyle=\sum_{i=1}^{m}\bm{X}^{\top}_{i}\dot{\bm{W}}_{i}\bm{X}_{i}-\hat{\bm{G}}^{-1}(\bm{X}^{\top}_{i}\dot{\bm{W}}_{i}\bm{X}_{i}+\hat{\bm{G}}^{-1})^{-1}\bm{X}^{\top}_{i}\dot{\bm{W}}_{i}\bm{X}_{i}$
$\displaystyle=\bm{B}_{1}-\sum_{i=1}^{m}\hat{\bm{G}}^{-1}(\bm{X}^{\top}_{i}\dot{\bm{W}}_{i}\bm{X}_{i}+\hat{\bm{G}}^{-1})^{-1}\bm{X}^{\top}_{i}\dot{\bm{W}}_{i}\bm{X}_{i}.$

It follows that

$\displaystyle\bm{C}$
$\displaystyle=\sum_{i=1}^{m}\hat{\bm{G}}^{-1}(\bm{X}^{\top}_{i}\dot{\bm{W}}_{i}\bm{X}_{i}+\hat{\bm{G}}^{-1})^{-1}\bm{X}^{\top}_{i}\dot{\bm{W}}_{i}\bm{X}_{i}=\sum_{i=1}^{m}\bm{X}^{\top}_{i}\dot{\bm{W}}_{i}\bm{X}_{i}(\bm{X}^{\top}_{i}\dot{\bm{W}}_{i}\bm{X}_{i}+\hat{\bm{G}}^{-1})^{-1}\hat{\bm{G}}^{-1},$
(S2.4)

where the second equality arises from the fact that as a covariance matrix,
$\bm{C}$ must be symmetric. We may also write $\bm{C}$ as

$\displaystyle\sum_{i=1}^{m}\hat{\bm{G}}^{-1}(\bm{X}^{\top}_{i}\dot{\bm{W}}_{i}\bm{X}_{i}+\hat{\bm{G}}^{-1})^{-1}\bm{X}^{\top}_{i}\dot{\bm{W}}_{i}\bm{X}_{i}$
$\displaystyle=\sum_{i=1}^{m}\\{\bm{I}_{p}-\hat{\bm{G}}^{-1}(\bm{X}^{\top}_{i}\dot{\bm{W}}_{i}\bm{X}_{i}+\hat{\bm{G}}^{-1})^{-1}\\}\hat{\bm{G}}^{-1}$
$\displaystyle=\hat{\bm{G}}^{-1}\sum_{i=1}^{m}\\{\bm{I}_{p}-(\bm{X}^{\top}_{i}\dot{\bm{W}}_{i}\bm{X}_{i}+\hat{\bm{G}}^{-1})^{-1}\hat{\bm{G}}^{-1}\\}.$
(S2.5)

Note that $\bm{C}$ is of order $O_{p}(m)$ component-wise in probability in
both the conditional and unconditional regimes. Using the fact that
$\bm{C}^{-1}$ must also be symmetric, we obtain

$\displaystyle\bm{C}^{-1}$
$\displaystyle=\left\\{\sum_{i=1}^{m}(\bm{X}^{\top}_{i}\dot{\bm{W}}_{i}\bm{X}_{i}+\hat{\bm{G}}^{-1})^{-1}\bm{X}^{\top}_{i}\dot{\bm{W}}_{i}\bm{X}_{i}\right\\}^{-1}\hat{\bm{G}}=\hat{\bm{G}}\left\\{\sum_{i=1}^{m}\bm{X}^{\top}_{i}\dot{\bm{W}}_{i}\bm{X}_{i}(\bm{X}^{\top}_{i}\dot{\bm{W}}_{i}\bm{X}_{i}+\hat{\bm{G}}^{-1})^{-1}\right\\}^{-1}$
(S2.6)

or equivalently

$\displaystyle\bm{C}^{-1}$
$\displaystyle=\bm{C}^{-1\top}=\left[\sum_{i=1}^{m}\\{\bm{I}_{p}-(\bm{X}^{\top}_{i}\dot{\bm{W}}_{i}\bm{X}_{i}+\hat{\bm{G}}^{-1})^{-1}\hat{\bm{G}}^{-1}\\}\right]^{-1}\hat{\bm{G}}$
$\displaystyle=\left\\{m^{-1}\bm{I}_{p}+m^{-1}\sum_{i=1}^{m}(\bm{X}^{\top}_{i}\dot{\bm{W}}_{i}\bm{X}_{i}+\hat{\bm{G}}^{-1})^{-1}\hat{\bm{G}}^{-1}\bm{C}^{-1}\right\\}\hat{\bm{G}},$
(S2.7)

where the last line is derived from (a special case of) the Woodbury identity,
given by
$(\bm{Q}-\bm{R})^{-1}=\bm{Q}^{-1}+\bm{Q}^{-1}\bm{R}(\bm{Q}-\bm{R})^{-1}$ for
arbitrary matrices $\bm{Q}$ and $\bm{R}$ such that $\bm{Q}$ and
$(\bm{Q}-\bm{R})$ are invertible. The first term in (S2.7) is the dominating
term, being of order $O(m^{-1})$, while the second term is
$O_{p}(m^{-1}n_{L}^{-1})$ in both the conditional and unconditional regimes.
We will use all the above forms of $\bm{C}$ and $\bm{C}^{-1}$ in subsequent
developments. Similarly, we can apply the Woodbury identity to
$(\bm{B}_{3}+\bm{B}_{4})^{-1}$ and
$(\bm{X}^{\top}_{i}\dot{\bm{W}}_{i}\bm{X}_{i}+\hat{\bm{G}}^{-1})^{-1}$ to
obtain
$n_{L}(\bm{B}_{3}+\bm{B}_{4})^{-1}=n_{L}\bm{B}_{3}^{-1}-n_{L}\bm{B}_{3}^{-1}\bm{B}_{4}(\bm{B}_{3}+\bm{B}_{4})^{-1}=O_{p}(1)+O_{p}(n_{L}^{-1})$
and
$n_{i}(\bm{X}^{\top}_{i}\dot{\bm{W}}_{i}\bm{X}_{i}+\hat{\bm{G}}^{-1})^{-1}=n_{i}(\bm{X}^{\top}_{i}\dot{\bm{W}}_{i}\bm{X}_{i})^{-1}-n_{i}(\bm{X}^{\top}_{i}\dot{\bm{W}}_{i}\bm{X}_{i})^{-1}\hat{\bm{G}}^{-1}(\bm{X}^{\top}_{i}\dot{\bm{W}}_{i}\bm{X}_{i}+\hat{\bm{G}}^{-1})^{-1}=O_{p}(1)+O_{p}(n_{i}^{-1})$,
where the order results hold component-wise. These hold irrespective of
whether we are conditioning on the random effects.

To further simplify expressions, for the rest of this article we will only use
order results when representing quantities associated with these smaller order
terms. Furthermore, as we want the derivations for the remainder of this
section to be applicable to both the conditional and unconditional regime, we
will not distinguish between $O(\cdot)$ and $O_{p}(\cdot)$ in the following
developments, and simply use $O_{p}()$ to represent both as appropriate. The
terms we use “big-O notation” for will have the same order under both the
conditional and unconditional regime. To simplify expressions, we will also
drop the dependence on $\bm{\theta}$, unless stated otherwise.

Finally, it is worth emphasising that

$\displaystyle[-\bm{I}_{p},\bm{I}_{p},\ldots,\bm{I}_{p}]\left[\begin{matrix}\dot{\phi}^{-1}\bm{X}^{\top}(\bm{y}-\dot{\bm{\mu}})\\\
\dot{\phi}^{-1}\bm{Z}^{\top}(\bm{y}-\dot{\bm{\mu}})\end{matrix}\right]=-\bm{S}_{1}+\sum_{i=1}^{m}\bm{S}_{2i}=\bm{S}_{1}-\sum_{i=1}^{m}\bm{S}_{2i}=\bm{0}_{p},$
(S2.8)

due to the $\bm{X}_{i}=\bm{Z}_{i}$ assumption. This is a key identity that is
critical to the proofs throughout this article.

We now use the expressions above to multiply out
$-\\{\nabla^{2}Q(\dot{\bm{\theta}})\\}^{-1}\nabla Q(\dot{\bm{\theta}})$ and
obtain expressions for $\hat{\bm{\beta}}-\dot{\bm{\beta}}$ and
$\hat{\bm{b}}-\dot{\bm{b}}$. From equation (S2.2), the first $p$ components of
$\hat{\bm{\theta}}-\dot{\bm{\theta}}$ are

$\displaystyle\hat{\bm{\beta}}-\dot{\bm{\beta}}$
$\displaystyle=\left[\begin{matrix}\bm{C}^{-1}&-\bm{C}^{-1}\bm{B}_{2}(\bm{B}_{3}+\bm{B}_{4})^{-1}\end{matrix}\right]\nabla
Q+\frac{1}{2}\\{\bm{B}^{-1}\bm{R}(\tilde{\bm{\theta}})\\}_{[1:p]}$
$\displaystyle=\bm{C}^{-1}\left[\begin{matrix}\bm{I}_{p}&-[\bm{I}_{p}-\hat{\bm{G}}^{-1}(\bm{X}^{\top}_{1}\dot{\bm{W}}_{1}\bm{X}_{1}+\hat{\bm{G}}^{-1})^{-1},\ldots,\bm{I}_{p}-\hat{\bm{G}}^{-1}(\bm{X}^{\top}_{m}\dot{\bm{W}}_{m}\bm{X}_{m}+\hat{\bm{G}}^{-1})^{-1}]\end{matrix}\right]\nabla
Q$
$\displaystyle+\frac{1}{2}\\{\bm{B}^{-1}\bm{R}(\tilde{\bm{\theta}})\\}_{[1:p]}$
$\displaystyle=\bm{C}^{-1}\left(\bm{S}_{1}-\sum_{i=1}^{m}\bm{S}_{2i}-\sum_{i=1}^{m}\bm{S}_{3i}\right)+\bm{C}^{-1}\hat{\bm{G}}^{-1}\Bigg{\\{}\sum_{i=1}^{m}(\bm{X}_{i}^{\top}\dot{\bm{W}}_{i}\bm{X}_{i}+\hat{\bm{G}}^{-1})^{-1}\bm{S}_{2i}$
$\displaystyle+\sum_{i=1}^{m}(\bm{X}_{i}^{\top}\dot{\bm{W}}_{i}\bm{X}_{i}+\hat{\bm{G}}^{-1})^{-1}\bm{S}_{3i}\Bigg{\\}}+\frac{1}{2}\\{\bm{B}^{-1}\bm{R}(\tilde{\bm{\theta}})\\}_{[1:p]}$
$\displaystyle=\bm{C}^{-1}\hat{\bm{G}}^{-1}\Bigg{\\{}\sum_{i=1}^{m}(\bm{X}_{i}^{\top}\dot{\bm{W}}_{i}\bm{X}_{i}+\hat{\bm{G}}^{-1})^{-1}\bm{S}_{2i}-\hat{\bm{G}}\sum_{i=1}^{m}\bm{S}_{3i}$
$\displaystyle+\sum_{i=1}^{m}(\bm{X}_{i}^{\top}\dot{\bm{W}}_{i}\bm{X}_{i}+\hat{\bm{G}}^{-1})^{-1}\bm{S}_{3i}\Bigg{\\}}+\frac{1}{2}\\{\bm{B}^{-1}\bm{R}(\tilde{\bm{\theta}})\\}_{[1:p]},$

where the final equality uses equation (S2.8). Thus, letting
$\bm{V}_{1}=\sum_{i=1}^{m}(\bm{X}_{i}^{\top}\dot{\bm{W}}_{i}\bm{X}_{i}+\hat{\bm{G}}^{-1})^{-1}\bm{S}_{2i}-\hat{\bm{G}}\sum_{i=1}^{m}\bm{S}_{3i}+\sum_{i=1}^{m}(\bm{X}_{i}^{\top}\dot{\bm{W}}_{i}\bm{X}_{i}+\hat{\bm{G}}^{-1})^{-1}\bm{S}_{3i}$
and applying equation (S2.7), we obtain

$\displaystyle\hat{\bm{\beta}}-\dot{\bm{\beta}}$
$\displaystyle=m^{-1}\bm{V}_{1}+\frac{1}{2}\\{\bm{B}^{-1}\bm{R}(\tilde{\bm{\theta}})\\}_{[1:p]}+O_{p}(n_{L}^{-1})\times
m^{-1}\bm{V}_{1}.$

Finally, using the Woodbury identity for
$(\bm{X}^{\top}_{i}\dot{\bm{W}}_{i}\bm{X}_{i}+\hat{\bm{G}}^{-1})^{-1}$, we
have that
$\sum_{i=1}^{m}(\bm{X}_{i}^{\top}\dot{\bm{W}}_{i}\bm{X}_{i}+\hat{\bm{G}}^{-1})^{-1}\bm{S}_{2i}=\sum_{i=1}^{m}(\bm{X}_{i}^{\top}\dot{\bm{W}}_{i}\bm{X}_{i})^{-1}\bm{S}_{2i}+\sum_{i=1}^{m}O_{p}(n_{L}^{-2})\bm{S}_{2i}$.
Letting
$\bm{V}_{2}=\sum_{i=1}^{m}(\bm{X}_{i}^{\top}\dot{\bm{W}}_{i}\bm{X}_{i})^{-1}\bm{S}_{2i}-\hat{\bm{G}}\sum_{i=1}^{m}\bm{S}_{3i}+\sum_{i=1}^{m}(\bm{X}_{i}^{\top}\dot{\bm{W}}_{i}\bm{X}_{i}+\hat{\bm{G}}^{-1})^{-1}\bm{S}_{3i}$,
we obtain

$\displaystyle\hat{\bm{\beta}}-\dot{\bm{\beta}}$
$\displaystyle=m^{-1}\bm{V}_{2}+\frac{1}{2}\\{\bm{B}^{-1}\bm{R}(\tilde{\bm{\theta}})\\}_{[1:p]}+O_{p}(n_{L}^{-1})\times
m^{-1}\bm{V}_{1}+m^{-1}\sum_{i=1}^{m}O_{p}(n_{L}^{-2})\bm{S}_{2i}.$

Next, the last $mp$ components of $\hat{\bm{\theta}}-\dot{\bm{\theta}}$ are

$\displaystyle\hat{\bm{b}}-\dot{\bm{b}}$
$\displaystyle=[-(\bm{B}_{3}+\bm{B}_{4})^{-1}\bm{B}_{2}^{\top}\bm{C}^{-1}\quad(\bm{B}_{3}+\bm{B}_{4})^{-1}+(\bm{B}_{3}+\bm{B}_{4})^{-1}\bm{B}_{2}^{\top}\bm{C}^{-1}\bm{B}_{2}(\bm{B}_{3}+\bm{B}_{4})^{-1}]\nabla
Q$
$\displaystyle+\frac{1}{2}\\{\bm{B}^{-1}\bm{R}(\tilde{\bm{\theta}})\\}_{[p+1:(m+1)p]}$
$\displaystyle=[-(\bm{B}_{3}+\bm{B}_{4})^{-1}\bm{B}_{2}^{\top}\bm{C}^{-1}\quad(\bm{B}_{3}+\bm{B}_{4})^{-1}\bm{B}_{2}^{\top}\bm{C}^{-1}\bm{B}_{2}(\bm{B}_{3}+\bm{B}_{4})^{-1}]\nabla
Q$ $\displaystyle+[\bm{0}_{mp\times p}\quad(\bm{B}_{3}+\bm{B}_{4})^{-1}]\nabla
Q+\frac{1}{2}\\{\bm{B}^{-1}\bm{R}(\tilde{\bm{\theta}})\\}_{[p+1:(m+1)p]}$
$\displaystyle=-(\bm{B}_{3}+\bm{B}_{4})^{-1}\bm{B}_{2}^{\top}[\bm{C}^{-1}\quad-\bm{C}^{-1}\bm{B}_{2}(\bm{B}_{3}+\bm{B}_{4})^{-1}]\nabla
Q$ $\displaystyle+[\bm{0}_{mp\times p}\quad(\bm{B}_{3}+\bm{B}_{4})^{-1}]\nabla
Q+\frac{1}{2}\\{\bm{B}^{-1}\bm{R}(\tilde{\bm{\theta}})\\}_{[p+1:(m+1)p]}.$

Notice that we already have an expression for
$[\bm{C}^{-1}\quad-\bm{C}^{-1}\bm{B}_{2}(\bm{B}_{3}+\bm{B}_{4})^{-1}]\nabla Q$
from the fixed effects above. Namely, it is
$m^{-1}\bm{V}_{1}+O_{p}(n_{L}^{-1})\times m^{-1}\bm{V}_{1}$. Thus we have

$\displaystyle\hat{\bm{b}}-\dot{\bm{b}}$
$\displaystyle=-(\bm{B}_{3}+\bm{B}_{4})^{-1}\bm{B}_{2}^{\top}(m^{-1}\bm{V}_{1}+O_{p}(n_{L}^{-1})\times
m^{-1}\bm{V}_{1})$
$\displaystyle+(\bm{B}_{3}+\bm{B}_{4})^{-1}\bm{S}_{6}+\frac{1}{2}\\{\bm{B}^{-1}\bm{R}(\tilde{\bm{\theta}})\\}_{[p+1:(m+1)p]}.$

Applying the Woodbury identity for $(\bm{B}_{3}+\bm{B}_{4})^{-1}$, we obtain

$\displaystyle\hat{\bm{b}}-\dot{\bm{b}}$
$\displaystyle=-\bm{1}_{m}\otimes(m^{-1}\bm{V}_{1}+O_{p}(n_{L}^{-1})\times
m^{-1}\bm{V}_{1})+O_{p}(n_{L}^{-1})(m^{-1}\bm{V}_{1}+O_{p}(n_{L}^{-1})\times
m^{-1}\bm{V}_{1})$
$\displaystyle+\bm{B}_{3}^{-1}\bm{S}_{6}+O_{p}(n_{L}^{-2})\bm{S}_{6}+\frac{1}{2}\\{\bm{B}^{-1}\bm{R}(\tilde{\bm{\theta}})\\}_{[p+1:(m+1)p]}$
$\displaystyle=-\bm{1}_{m}\otimes m^{-1}\bm{V}_{1}+O_{p}(n_{L}^{-1})\times
m^{-1}\bm{V}_{1}+O_{p}(n_{L}^{-2})\times m^{-1}\bm{V}_{1}$
$\displaystyle+\bm{B}_{3}^{-1}\bm{S}_{4}+\bm{B}_{3}^{-1}\bm{S}_{5}+O_{p}(n_{L}^{-2})\bm{S}_{6}+\frac{1}{2}\\{\bm{B}^{-1}\bm{R}(\tilde{\bm{\theta}})\\}_{[p+1:(m+1)p]}.$

Replacing all the $\bm{V_{\cdot}}$ and $\bm{S_{\cdot}}$ terms in the above
with their definitions, we finally obtain

$\displaystyle\hat{\bm{\beta}}-\dot{\bm{\beta}}$
$\displaystyle=m^{-1}\sum_{i=1}^{m}(\bm{X}_{i}^{\top}\dot{\bm{W}}_{i}\bm{X}_{i})^{-1}\dot{\phi}^{-1}\bm{X}_{i}^{\top}(\bm{y}_{i}-\dot{\bm{\mu}}_{i})+m^{-1}\sum_{i=1}^{m}\dot{\bm{b}}_{i}$
$\displaystyle-m^{-1}\sum_{i=1}^{m}(\bm{X}_{i}^{\top}\dot{\bm{W}}_{i}\bm{X}_{i}+\hat{\bm{G}}^{-1})^{-1}\hat{\bm{G}}^{-1}\dot{\bm{b}}_{i}+\frac{1}{2}\\{\bm{B}^{-1}\bm{R}(\tilde{\bm{\theta}})\\}_{[1:p]}$
$\displaystyle+O_{p}(n_{L}^{-1})\Bigg{\\{}m^{-1}\sum_{i=1}^{m}(\bm{X}_{i}^{\top}\dot{\bm{W}}_{i}\bm{X}_{i}+\hat{\bm{G}}^{-1})^{-1}\dot{\phi}^{-1}\bm{X}_{i}^{\top}(\bm{y}_{i}-\dot{\bm{\mu}}_{i})+m^{-1}\sum_{i=1}^{m}\dot{\bm{b}}_{i}$
$\displaystyle-m^{-1}\sum_{i=1}^{m}(\bm{X}_{i}^{\top}\dot{\bm{W}}_{i}\bm{X}_{i}+\hat{\bm{G}}^{-1})^{-1}\hat{\bm{G}}^{-1}\dot{\bm{b}}_{i}\Bigg{\\}}+m^{-1}\sum_{i=1}^{m}O_{p}(n_{L}^{-2})\dot{\phi}^{-1}\bm{X}_{i}^{\top}(\bm{y}_{i}-\dot{\bm{\mu}}_{i}),$
(S2.9)

and

$\displaystyle\hat{\bm{b}}-\dot{\bm{b}}$
$\displaystyle=-\bm{1}_{m}\otimes\Bigg{\\{}m^{-1}\sum_{i=1}^{m}(\bm{X}_{i}^{\top}\dot{\bm{W}}_{i}\bm{X}_{i}+\hat{\bm{G}}^{-1})^{-1}\dot{\phi}^{-1}\bm{X}_{i}^{\top}(\bm{y}_{i}-\dot{\bm{\mu}}_{i})$
$\displaystyle+m^{-1}\sum_{i=1}^{m}\dot{\bm{b}}_{i}-m^{-1}\sum_{i=1}^{m}(\bm{X}_{i}^{\top}\dot{\bm{W}}_{i}\bm{X}_{i}+\hat{\bm{G}}^{-1})^{-1}\hat{\bm{G}}^{-1}\dot{\bm{b}}_{i}\Bigg{\\}}$
$\displaystyle+\bm{B}_{3}^{-1}\\{\dot{\phi}^{-1}\bm{Z}^{\top}(\bm{y}-\dot{\bm{\mu}})\\}-\bm{B}_{3}^{-1}\\{(\bm{I}_{m}\otimes\hat{\bm{G}}^{-1})\dot{\bm{b}}\\}+\frac{1}{2}\\{\bm{B}^{-1}\bm{R}(\tilde{\bm{\theta}})\\}_{[p+1:(m+1)p]}$
$\displaystyle+O_{p}(n_{L}^{-1})\Bigg{\\{}m^{-1}\sum_{i=1}^{m}(\bm{X}_{i}^{\top}\dot{\bm{W}}_{i}\bm{X}_{i}+\hat{\bm{G}}^{-1})^{-1}\dot{\phi}^{-1}\bm{X}_{i}^{\top}(\bm{y}_{i}-\dot{\bm{\mu}}_{i})$
$\displaystyle+m^{-1}\sum_{i=1}^{m}\dot{\bm{b}}_{i}-m^{-1}\sum_{i=1}^{m}(\bm{X}_{i}^{\top}\dot{\bm{W}}_{i}\bm{X}_{i}+\hat{\bm{G}}^{-1})^{-1}\hat{\bm{G}}^{-1}\dot{\bm{b}}_{i}\Bigg{\\}}$
$\displaystyle+O_{p}(n_{L}^{-2})\Bigg{\\{}m^{-1}\sum_{i=1}^{m}(\bm{X}_{i}^{\top}\dot{\bm{W}}_{i}\bm{X}_{i}+\hat{\bm{G}}^{-1})^{-1}\dot{\phi}^{-1}\bm{X}_{i}^{\top}(\bm{y}_{i}-\dot{\bm{\mu}}_{i})$
$\displaystyle+m^{-1}\sum_{i=1}^{m}\dot{\bm{b}}_{i}-m^{-1}\sum_{i=1}^{m}(\bm{X}_{i}^{\top}\dot{\bm{W}}_{i}\bm{X}_{i}+\hat{\bm{G}}^{-1})^{-1}\hat{\bm{G}}^{-1}\dot{\bm{b}}_{i}\Bigg{\\}}$
$\displaystyle+O_{p}(n_{L}^{-2})\\{\dot{\phi}^{-1}\bm{Z}^{\top}(\bm{y}-\dot{\bm{\mu}})-(\bm{I}_{m}\otimes\hat{\bm{G}}^{-1})\dot{\bm{b}}\\}.$
(S2.10)

The expressions for $\hat{\bm{\beta}}-\dot{\bm{\beta}}$ and
$\hat{\bm{b}}-\dot{\bm{b}}$ above underlie our proofs. We use these same
expressions in both the conditional and unconditional regimes, but the
asymptotic behaviours of the terms on the right hand side, and the way we
treat them, will differ greatly between the two cases.

As we will show later, the key leading terms for the fixed effects are
$m^{-1}\sum_{i=1}^{m}(\bm{X}_{i}^{\top}\dot{\bm{W}}_{i}\bm{X}_{i})^{-1}\dot{\phi}^{-1}\bm{X}_{i}^{\top}(\bm{y}_{i}-\dot{\bm{\mu}}_{i})$
and $m^{-1}\sum_{i=1}^{m}\dot{\bm{b}}_{i}$. The key leading terms for the
random effects are $-\bm{1}_{m}\otimes m^{-1}\sum_{i=1}^{m}\dot{\bm{b}}_{i}$
and
$\bm{B}_{3}^{-1}\\{\dot{\phi}^{-1}\bm{Z}^{\top}(\bm{y}-\dot{\bm{\mu}})\\}$.
When conditioning on the random effects $\dot{\bm{b}}$, we have
$m^{-1}\sum_{i=1}^{m}\dot{\bm{b}}_{i}=O(1)$, while in the unconditional regime
the same quantity is of order $O_{p}(m^{-1/2})$ in probability. In both the
conditional and unconditional regimes, we have that
$m^{-1}\sum_{i=1}^{m}(\bm{X}_{i}^{\top}\dot{\bm{W}}_{i}\bm{X}_{i})^{-1}\dot{\phi}^{-1}\bm{X}_{i}^{\top}(\bm{y}_{i}-\dot{\bm{\mu}}_{i})$
is of order $O_{p}(N^{-1/2})$ component-wise, while the quantity
$\bm{B}_{3}^{-1}\\{\dot{\phi}^{-1}\bm{Z}^{\top}(\bm{y}-\dot{\bm{\mu}})\\}$ is
of order $O_{p}(n_{L}^{-1/2})$ component-wise.

###  S2.1 Proof of Theorem 1

The dominating terms on the right hand sides of equations (S2.9) and (S2.10)
are $m^{-1}\sum_{i=1}^{m}\dot{\bm{b}}_{i}$ and $\bm{1}_{m}\otimes
m^{-1}\sum_{i=1}^{m}\dot{\bm{b}}_{i}$ for the fixed and random effects,
respectively. Conditional on the random effects $\dot{\bm{b}}_{i}$, these
dominating terms are deterministic and of order $O(1)$. Thus we treat them as
bias terms and move them to the left hand side. Next, by Conditions (C1)-(C2),
$\bm{B}_{3}^{-1}$ is a component-wise $O(n_{L}^{-1})$ block-diagonal matrix,
while we also have $\bm{B}_{2}^{\top}=O(n_{U})$,
$\bm{X}_{i}^{\top}\dot{\bm{W}}_{i}\bm{X}_{i}^{\top}=O(n_{i})$, and
$\bm{C}^{-1}=O(m^{-1})$ component-wise. Since
$E\\{\bm{Z}^{\top}(\bm{y}-\dot{\bm{\mu}})|\dot{\bm{b}}\\}=\bm{0}_{mp}$ and
$\text{Var}\\{\bm{Z}^{\top}(\bm{y}-\dot{\bm{\mu}})|\dot{\bm{b}}\\}=\bm{Z}^{\top}\dot{\bm{W}}\bm{Z}$,
we obtain
$\dot{\phi}^{-1}\bm{D}_{r}^{-1}\bm{Z}^{\top}(\bm{y}-\dot{\bm{\mu}})=O_{p}(1)$
using Chebyshev’s inequality and the conditional independence.

Multiplying both sides of (S2.9) and (S2.10) by $N^{1/2}$ and $\bm{D}_{r}$
respectively, and applying the order results for the remainder term in Section
S3.1, we obtain

$\displaystyle
N^{1/2}\left(\hat{\bm{\beta}}-\dot{\bm{\beta}}-m^{-1}\sum_{i=1}^{m}\dot{\bm{b}}_{i}\right)$
$\displaystyle=m^{-1/2}\sum_{i=1}^{m}n^{1/2}n_{i}^{-1/2}(n_{i}^{-1}\bm{X}_{i}^{\top}\dot{\bm{W}}_{i}\bm{X}_{i})^{-1}n_{i}^{-1/2}\dot{\phi}^{-1}\bm{X}_{i}^{\top}(\bm{y}_{i}-\dot{\bm{\mu}}_{i})$
$\displaystyle+O_{p}(m^{1/2}n_{L}^{-1/2}),$

and

$\displaystyle\bm{D}_{r}\left(\hat{\bm{b}}-\dot{\bm{b}}+\bm{1}_{m}\otimes
m^{-1}\sum_{i=1}^{m}\dot{\bm{b}}_{i}\right)$
$\displaystyle=\bm{D}_{r}\bm{B}_{3}^{-1}\bm{D}_{r}\bm{D}_{r}^{-1}\\{\dot{\phi}^{-1}\bm{Z}^{\top}(\bm{y}-\dot{\bm{\mu}})\\}+O_{p}(n_{L}^{-1/2}).$

Recalling that $\bm{X}_{i}=\bm{Z}_{i}$, to prove Theorem 1 we will show a
Lindeberg condition for

$\displaystyle\bm{A}\left[\begin{matrix}m^{-1/2}\sum_{i=1}^{m}n^{1/2}n_{i}^{-1/2}(n_{i}^{-1}\bm{X}_{i}^{\top}\dot{\bm{W}}_{i}\bm{X}_{i})^{-1}n_{i}^{-1/2}\dot{\phi}^{-1}\bm{X}_{i}^{\top}(\bm{y}_{i}-\dot{\bm{\mu}}_{i})\\\
(n_{1}^{-1}\bm{X}_{1}^{\top}\dot{\bm{W}}_{1}\bm{X}_{1})^{-1}\\{n_{1}^{-1/2}\dot{\phi}^{-1}\bm{X}_{1}^{\top}(\bm{y}_{1}-\dot{\bm{\mu}}_{1})\\}\\\
\vdots\\\
(n_{m}^{-1}\bm{X}_{m}^{\top}\dot{\bm{W}}_{m}\bm{X}_{m})^{-1}\\{n_{m}^{-1/2}\dot{\phi}^{-1}\bm{X}_{m}^{\top}(\bm{y}_{m}-\dot{\bm{\mu}}_{m})\\}\end{matrix}\right]=:S,$

and thus apply the Lindeberg-Feller central limit theorem, from which the
result follows from Slutsky’s theorem.

To prove the condition, first define
$\bm{U}=[\bm{Z}\bm{B}_{3}^{-1}(\bm{1}_{m}\otimes\bm{I}_{p}),\bm{Z}\bm{B}_{3}^{-1}]$,
and $\bm{U}_{k}$ as the $k$th row of $\bm{U}$, noting it only has $2p$ non-
zero components. Then we can write
$S=\sum_{k=1}^{N}\bm{A}\bm{D}\bm{U}_{k}\dot{\phi}^{-1}\\{y_{k}-\mu_{k}(\dot{\bm{\theta}})\\}\overset{\Delta}{=}\sum_{k=1}^{N}\bm{\xi}_{k}$,
where $y_{k}$ is the $k$th component in
$(y_{11},y_{12},\ldots,y_{1n_{1}},y_{21},\ldots,y_{mn_{m}})^{\top}$, and
similarly for $\mu_{k}(\dot{\bm{\theta}})$.

Conditional on $\dot{\bm{b}}$, the quantities $\\{\bm{\xi}_{k}\\}_{k=1}^{N}$
are independent $q$-vectors with expectation zero and covariance
$\text{Var}(\bm{\xi}_{k}|\dot{\bm{b}})=\bm{A}\bm{D}\bm{U}_{k}W_{k}\bm{U}_{k}^{\top}\bm{D}\bm{A}^{\top}$,
where $W_{k}$ is the $k$th diagonal component in $\dot{\bm{W}}$. Therefore, we
have that

$\displaystyle\sum_{k=1}^{N}\text{Var}(\bm{\xi}_{k}|\dot{\bm{b}})=\sum_{k=1}^{N}\bm{A}\bm{D}\bm{U}_{k}W_{k}\bm{U}_{k}^{\top}\bm{D}\bm{A}^{\top}$
$\displaystyle=\bm{A}\left[\begin{matrix}\frac{1}{m}\sum_{i=1}^{m}\frac{n}{n_{i}}\left(\frac{\bm{X}^{\top}_{i}\dot{\bm{W}}_{i}\bm{X}_{i}}{n_{i}}\right)^{-1}&\frac{1}{\sqrt{m}}\sqrt{\frac{n}{n_{1}}}\left(\frac{\bm{X}^{\top}_{1}\dot{\bm{W}}_{1}\bm{X}_{1}}{n_{1}}\right)^{-1}&\cdots&\frac{1}{\sqrt{m}}\sqrt{\frac{n}{n_{m}}}\left(\frac{\bm{X}^{\top}_{m}\dot{\bm{W}}_{m}\bm{X}_{m}}{n_{m}}\right)^{-1}\\\
\frac{1}{\sqrt{m}}\sqrt{\frac{n}{n_{1}}}\left(\frac{\bm{X}^{\top}_{1}\dot{\bm{W}}_{1}\bm{X}_{1}}{n_{1}}\right)^{-1}&\left(\frac{\bm{X}^{\top}_{1}\dot{\bm{W}}_{1}\bm{X}_{1}}{n_{1}}\right)^{-1}&\bm{0}&\bm{0}\\\
\vdots&\bm{0}&\ddots&\bm{0}\\\
\frac{1}{\sqrt{m}}\sqrt{\frac{n}{n_{m}}}\left(\frac{\bm{X}^{\top}_{m}\dot{\bm{W}}_{m}\bm{X}_{m}}{n_{m}}\right)^{-1}&\bm{0}&\bm{0}&\left(\frac{\bm{X}^{\top}_{m}\dot{\bm{W}}_{m}\bm{X}_{m}}{n_{m}}\right)^{-1}\end{matrix}\right]\bm{A}^{\top}.$

Hence using the finite selection property of $\bm{A}$, and the fact that
${m}^{-1/2}n^{1/2}n_{i}^{-1/2}\left(n_{i}^{-1}\bm{X}^{\top}_{i}\dot{\bm{W}}_{i}\bm{X}_{i}\right)^{-1}=o(1)$
component-wise, we obtain

$\displaystyle\underset{m,n_{L}\rightarrow\infty}{\lim}\sum_{k=1}^{N}\text{Cov}(\bm{\xi}_{k}|\dot{\bm{b}})$
$\displaystyle=\underset{m,n_{L}\rightarrow\infty}{\lim}\bm{A}\,\text{bdiag}\Bigg{\\{}\frac{1}{m}\sum_{i=1}^{m}\frac{n}{n_{i}}\left(\frac{\bm{X}^{\top}_{i}\dot{\bm{W}}_{i}\bm{X}_{i}}{n_{i}}\right)^{-1},\left(\frac{\bm{X}^{\top}_{1}\dot{\bm{W}}_{1}\bm{X}_{1}}{n_{1}}\right)^{-1},\ldots,\left(\frac{\bm{X}^{\top}_{m}\dot{\bm{W}}_{m}\bm{X}_{m}}{n_{m}}\right)^{-1}\Bigg{\\}}\bm{A}^{\top}$
$\displaystyle=\bm{\Omega}.$

Next, by the Cauchy-Schwarz inequality, we have

$\displaystyle
E\\{\|\bm{\xi}_{k}\|^{2}I(\|\bm{\xi}_{k}\|>\epsilon)|\dot{\bm{b}}\\}\leq
E(\|\bm{\xi}_{k}\|^{4}|\dot{\bm{b}})^{1/2}P(\|\bm{\xi}_{k}\|>\epsilon|\dot{\bm{b}})^{1/2}.$

Finally, we make a note about the form of
Cov$[\bm{D}\bm{U}_{k}\\{y_{k}-\mu_{k}(\dot{\bm{\theta}})\\}]$. Without loss of
generality, suppose $k=1$. Then

$\displaystyle\text{Cov}[\bm{D}\bm{U}_{1}\\{y_{1}-\mu_{1}(\dot{\bm{\theta}})\\}]=$
$\displaystyle\left[\begin{matrix}n(mn_{1}^{2})^{-1}\bm{H}_{1}\bm{x}_{11}W_{1}\bm{x}_{11}^{\top}\bm{H}_{1}^{\top}&n_{1}^{-1}m^{-1/2}(nn_{1}^{-1})^{1/2}\bm{H}_{1}\bm{x}_{11}W_{1}\bm{x}_{11}^{\top}\bm{H}_{1}^{\top}&\bm{0}\\\
n_{1}^{-1}m^{-1/2}(nn_{1}^{-1})^{1/2}\bm{H}_{1}\bm{x}_{11}W_{1}\bm{x}_{11}^{\top}\bm{H}_{1}^{\top}&n_{1}^{-1}\bm{H}_{1}\bm{x}_{11}W_{1}\bm{x}_{11}^{\top}\bm{H}_{1}^{\top}&\bm{0}\\\
\bm{0}&\bm{0}&\bm{0}\\\ \end{matrix}\right].$ (S2.11)

Again without loss of generality, consider the case
$\bm{A}=[\bm{I}_{2p},\bm{0}_{(p+p)\times(m-1)p}]$. Then by equation (S2.11)
and Chebyshev’s inequality, when $k\in\\{1,2,\ldots,n_{1}\\}$ we have that
$P(\|\bm{\xi}_{k}\|>\epsilon|\dot{\bm{b}})\leq\text{tr}\\{\text{Cov}(\bm{\xi}_{k}|\dot{\bm{b}})\\}/\epsilon^{2}=O(n_{1}^{-1})$.
Thus, given $\dot{\bm{b}}$, we obtain $\|\bm{\xi}_{k}\|=O_{p}(n_{1}^{-1/2})$
and $E(\|\bm{\xi}_{k}\|^{4}|\dot{\bm{b}})=O(n_{1}^{-2})$ by Conditions
(C1)-(C3) and the properties of the exponential family. However when
$k>n_{1}$, by equation (S2.11) and Chebyshev’s inequality, we have that
$P(\|\bm{\xi}_{k}\|>\epsilon|\dot{\bm{b}})\leq\text{tr}\\{\text{Cov}(\bm{\xi}_{k}|\dot{\bm{b}})\\}/\epsilon^{2}=O(N^{-1})$
since $n(mn_{1}^{2})^{-1}=O(N^{-1})$. Thus given $\dot{\bm{b}}$, it holds that
$\|\bm{\xi}_{k}\|=O_{p}(N^{-1/2})$ and
$E(\|\bm{\xi}_{k}\|^{4}|\dot{\bm{b}})=O(N^{-2})$. Therefore

$\displaystyle\sum_{k=1}^{N}E\\{\|\bm{\xi}_{k}\|^{2}I(\|\bm{\xi}_{k}\|>\epsilon)|\dot{\bm{b}}\\}$
$\displaystyle\leq\sum_{k=1}^{N}E(\|\bm{\xi}_{k}\|^{4}|\dot{\bm{b}})^{1/2}P(\|\bm{\xi}_{k}\|>\epsilon|\dot{\bm{b}})^{1/2}$
$\displaystyle=\sum_{k=1}^{n_{1}}E(\|\bm{\xi}_{k}\|^{4}|\dot{\bm{b}})^{1/2}P(\|\bm{\xi}_{k}\|>\epsilon|\dot{\bm{b}})^{1/2}$
$\displaystyle+\sum_{k=n_{1}+1}^{N}E(\|\bm{\xi}_{k}\|^{4}|\dot{\bm{b}})^{1/2}P(\|\bm{\xi}_{k}\|>\epsilon|\dot{\bm{b}})^{1/2}$
$\displaystyle\leq n_{1}\,\underset{1\leq k\leq
n_{1}}{\max}\\{E(\|\bm{\xi}_{k}\|^{4}|\dot{\bm{b}})^{1/2}P(\|\bm{\xi}_{k}\|>\epsilon|\dot{\bm{b}})^{1/2}\\}$
$\displaystyle+(N-n_{1})\,\underset{k>n_{1}}{\sup}\\{E(\|\bm{\xi}_{k}\|^{4}|\dot{\bm{b}})^{1/2}P(\|\bm{\xi}_{k}\|>\epsilon|\dot{\bm{b}})^{1/2}\\}$
$\displaystyle=n_{1}\times O(n_{1}^{-3/2})+(N-n_{1})\times O(N^{-3/2})$
$\displaystyle=O(n_{1}^{-1/2})+O(N^{-1/2})$ $\displaystyle=o(1).$

The required result follows by Conditions (C1)-(C2) and the Lindeberg-Feller
Central Limit Theorem. Furthermore, the general result holds straightforwardly
by replacing $n_{1}$ with $O(n_{L})$ in the above argument, noting that any
row of $\bm{A}$ can only select a fixed number of clusters.

###  S2.2 Proof of Equation (4)

For the Poisson pure random intercept model, we have
$\bm{B}=\text{diag}(ne^{\dot{b}_{1}}+\hat{\sigma}_{b}^{-2},\ldots,ne^{\dot{b}_{m}}+\hat{\sigma}_{b}^{-2})$
and
$\bm{R}(\tilde{\bm{\theta}})=\\{ne^{\tilde{b}_{1}}(\hat{b}_{1}-\dot{b}_{1})^{2},\ldots,ne^{\tilde{b}_{m}}(\hat{b}_{m}-\dot{b}_{m})^{2}\\}^{\top}$.
Next, suppose that $\bm{A}$ picks out the first random intercept, i.e.,
$\bm{A}=[1,\bm{0}_{m-1}^{\top}]$. Then we have

$\displaystyle{n}^{1/2}(\hat{b}_{1}-\dot{b}_{1})$
$\displaystyle={n}^{1/2}\bm{AB}^{-1}\nabla
Q(\dot{\bm{\theta}})+\frac{1}{2}{n}^{1/2}\bm{A}\bm{B}^{-1}\bm{R}(\tilde{\bm{\theta}})$
$\displaystyle=n^{-1/2}\left\\{\left(\sum_{j=1}^{n}y_{1j}-e^{\dot{b}_{1}}\right)-\dot{b}_{1}/\hat{\sigma}_{b}^{2}\right\\}\Bigg{/}\Big{\\{}e^{\dot{b}_{1}}+1/(\hat{\sigma}_{b}^{2}n)\Big{\\}}$
$\displaystyle-\frac{1}{2}\left\\{n^{1/2}e^{\tilde{b}_{1}}(\hat{b}_{1}-\dot{b}_{1})^{2}\right\\}\Bigg{/}\Big{\\{}e^{\dot{b}_{1}}+1/(\hat{\sigma}_{b}^{2}n)\Big{\\}}$
$\displaystyle=\Bigg{[}n^{-1/2}\left\\{\left(\sum_{j=1}^{n}y_{1j}-e^{\dot{b}_{1}}\right)-\dot{b}_{1}/\hat{\sigma}_{b}^{2}\right\\}\Bigg{/}\Big{\\{}e^{\dot{b}_{1}}+1/(\hat{\sigma}_{b}^{2}n)\Big{\\}}\Bigg{]}\Bigg{/}$
$\displaystyle\left[1+\left\\{\frac{1}{2}e^{\tilde{b}_{1}}(\hat{b}_{1}-\dot{b}_{1})\right\\}\Bigg{/}\Big{\\{}e^{\dot{b}_{1}}+1/(\hat{\sigma}_{b}^{2}n)\Big{\\}}\right]$
$\displaystyle=n^{-1/2}\sum_{j=1}^{n}(y_{1j}e^{-\dot{b}_{1}}-1)+o_{p}(1),$

where $\tilde{b}_{1}$ lies between $\hat{b}_{1}$ and $\dot{b}_{1}$, and for
the last line we have used the fact that $\hat{b}_{1}-\dot{b}_{1}=o_{p}(1)$.

Now, $\\{y_{1j}e^{-\dot{b}_{1}}-1\\}_{j=1}^{m}$ is an exchangeable collection
of uncorrelated random variables with mean zero and finite non-zero variance.
Furthermore, we have for $k\neq l$

$\displaystyle\text{Cov}\\{(y_{1k}e^{-\dot{b}_{1}}-1)^{2},(y_{1l}e^{-\dot{b}_{1}}-1)^{2}\\}$
$\displaystyle=E[\text{Cov}\\{(y_{1k}e^{-\dot{b}_{1}}-1)^{2},(y_{1l}e^{-\dot{b}_{1}}-1)^{2}|\dot{b}_{1}\\}]$
$\displaystyle\quad+\text{Cov}[E\\{(y_{1k}e^{-\dot{b}_{1}}-1)^{2}|\dot{b}_{1}\\},E\\{(y_{1l}e^{-\dot{b}_{1}}-1)^{2}|\dot{b}_{1}\\}]$
$\displaystyle=0+\text{Cov}(e^{-\dot{b}_{1}},e^{-\dot{b}_{1}})$
$\displaystyle=e^{\dot{\sigma}_{b}^{2}}(e^{\dot{\sigma}_{b}^{2}}-1)\neq 0.$

Thus by the Central Limit Theorem for exchangeable random variables [33], it
holds that
$n^{-1/2}\sum_{j=1}^{n}(y_{1j}e^{-\dot{b}_{1}}-1)\overset{D}{\not\rightarrow}N(0,e^{\dot{\sigma}_{b}^{2}})$.
Since we know
$\text{Var}\\{n^{-1/2}\sum_{j=1}^{n}(y_{1j}e^{-\dot{b}_{1}}-1)\\}=e^{\dot{\sigma}_{b}^{2}/2}$
and also that $n^{-1/2}\sum_{j=1}^{n}(y_{1j}e^{-\dot{b}_{1}}-1)=O_{p}(1)$ by
Chebyshev’s inequality, there is no other normalization possible for an
asymptotic normality result to hold.

Finally, we also have

$\displaystyle{n}^{1/2}(\hat{b}_{1}-\dot{b}_{1})$
$\displaystyle=n^{-1/2}\sum_{j=1}^{n}(y_{1j}e^{-\dot{b}_{1}}-1)+O_{p}(n^{-1/2})$
$\displaystyle\implies\hat{b}_{1}$
$\displaystyle=\dot{b}_{1}+n^{-1}\sum_{j=1}^{n}(y_{1j}e^{-\dot{b}_{1}}-1)+O_{p}(n^{-1})$
$\displaystyle=\dot{b}_{1}+o_{p}(1),\quad\text{by the Weak Law of Large
Numbers}.$

###  S2.3 Proof of Theorem 2

We begin by developing two key equations, (S2.12) and (S2.13), that will be
used throughout the unconditional regime. These are derived from equations
(S2.9) and (S2.10) and are used in the proofs of Theorems 2-5 as well as
Corollary 1. Under Conditions (C1)-(C2), the following order results are used:
$\bm{B}_{3}^{-1}$ is a component-wise $O_{p}(n_{L}^{-1})$ block-diagonal
matrix, $\bm{B}_{2}=O_{p}(n_{U})$ component-wise,
$\bm{X}_{i}^{\top}\dot{\bm{W}}_{i}\bm{X}_{i}^{\top}=O_{p}(n_{i})$ component-
wise, and $\bm{C}^{-1}=O_{p}(m^{-1})$ component-wise. Also, by the conditional
independence, we have

$\displaystyle
E\\{\bm{Z}^{\top}(\bm{y}-\dot{\bm{\mu}})\\}=E[E\\{\bm{Z}^{\top}(\bm{y}-\dot{\bm{\mu}})|\dot{\bm{b}}\\}]=\bm{0}_{mp},$
$\displaystyle\text{Var}\\{\bm{Z}^{\top}(\bm{y}-\dot{\bm{\mu}})\\}=E[\text{Var}\\{\bm{Z}^{\top}(\bm{y}-\dot{\bm{\mu}})|\dot{\bm{b}}\\}]+\text{Var}[E\\{\bm{Z}^{\top}(\bm{y}-\dot{\bm{\mu}})|\dot{\bm{b}}\\}]=E(\bm{Z}^{\top}\dot{\bm{W}}\bm{Z}),$

so that
$\dot{\phi}^{-1}\bm{D}_{r}^{-1}\bm{Z}^{\top}(\bm{y}-\dot{\bm{\mu}})=O_{p}(1)$
using Chebyshev’s inequality. Therefore we have the key equations

$\displaystyle\hat{\bm{\beta}}-\dot{\bm{\beta}}$
$\displaystyle=m^{-1}\sum_{i=1}^{m}\dot{\bm{b}}_{i}+O_{p}(N^{-1/2})+O_{p}(n_{L}^{-1})+\frac{1}{2}\\{\bm{B}^{-1}\bm{R}(\tilde{\bm{\theta}})\\}_{[1:p]}$
(S2.12)

and

$\displaystyle\hat{\bm{b}}-\dot{\bm{b}}$ $\displaystyle=-\bm{1}_{m}\otimes
m^{-1}\sum_{i=1}^{m}\dot{\bm{b}}_{i}+\bm{B}_{3}^{-1}\\{\dot{\phi}^{-1}\bm{Z}^{\top}(\bm{y}-\dot{\bm{\mu}})\\}$
$\displaystyle+O_{p}(N^{-1/2})+O_{p}(n_{L}^{-1})+\frac{1}{2}\\{\bm{B}^{-1}\bm{R}(\tilde{\bm{\theta}})\\}_{[p+1:(m+1)p]}.$
(S2.13)

By equation (S2.12), we have

$\displaystyle m^{1/2}(\hat{\bm{\beta}}-\dot{\bm{\beta}})$
$\displaystyle=m^{-1/2}\sum_{i=1}^{m}\dot{\bm{b}}_{i}+O_{p}(n_{L}^{-1/2})+O_{p}(m^{1/2}n_{L}^{-1})+\frac{1}{2}m^{1/2}\\{\bm{B}^{-1}\bm{R}(\tilde{\bm{\theta}})\\}_{[1:p]}.$

Next, we consider two separate scenarios. First, suppose that
$mn_{U}^{-1}\rightarrow\infty$. Then by the order results for the remainder
term in Section S3.2, the first $p$ components of
$\bm{D}^{*}\bm{B}^{-1}\bm{R}(\tilde{\bm{\theta}})$ are of order
$O_{p}(m^{1/2}n_{L}^{-1})$, and so the first $p$ components of
$\bm{D}^{*}(\hat{\bm{\theta}}-\dot{\bm{\theta}})$ can be shown to be

$\displaystyle{m}^{1/2}(\hat{\bm{\beta}}-\dot{\bm{\beta}})$
$\displaystyle=m^{-1/2}\sum_{i=1}^{m}\dot{\bm{b}}_{i}+o_{p}(1).$

The required result then follows from the independence of the random effects
and the normal assumption on the $\dot{\bm{b}}_{i}$; note the
$mn_{L}^{-2}\rightarrow 0$ assumption is required for the remainder term to be
smaller order than the linear term.

On the other hand, when $mn_{L}^{-1}\rightarrow 0$, the only difference from
the $mn_{L}^{-1}\rightarrow\infty$ case is that the first $p$ components of
$\bm{D}^{+}\bm{B}^{-1}\bm{R}(\tilde{\bm{\theta}})$ are now of order
$O_{p}(m^{-1/2})$ due to the different convergence rate of the prediction gap.
The result however follows along similar lines as above.

###  S2.4 Proof of Theorem 3

Again we consider two different scenarios. First, suppose
$mn_{L}^{-1}\rightarrow\infty$. Then from equation (S2.13) and the order
results for the remainder term in Section S3.2, we have that

$\displaystyle\bm{D}_{r}(\hat{\bm{b}}-\dot{\bm{b}})$
$\displaystyle=O_{p}(n_{U}^{1/2}m^{-1/2})+O_{p}(1)+O_{p}(n_{L}^{-1/2}).$

Based on the above, we obtain
$\bm{D}_{r}\hat{\bm{b}}=\bm{D}_{r}\dot{\bm{b}}+O_{p}(1)$, and thus
$\hat{\bm{b}}=\dot{\bm{b}}+O_{p}(n_{L}^{-1/2})$. The required result follows
by multiplying both sides by $\bm{A}_{r}$.

On the other hand, suppose now $mn_{L}^{-1}\rightarrow 0$. Then a
normalization by ${m}^{1/2}$ is needed instead, and the third derivative term
is consequently of order $O_{p}(m^{-1/2})$ in probability. We thus obtain

$\displaystyle{m}^{1/2}(\hat{\bm{b}}-\dot{\bm{b}})$
$\displaystyle=O_{p}(1)+O_{p}(m^{1/2}n_{L}^{-1/2})+O_{p}(m^{-1/2}),$

and the result follows.

As an side remark, note from the above proof when $mn_{L}^{-1}\rightarrow 0$,
it holds that $\|\hat{\bm{b}}-\dot{\bm{b}}\|_{2}=O_{p}(1)$, where
$\|\cdot\|_{2}$ denotes the $l_{2}$-norm. But if
$mn_{U}^{-1}\rightarrow\infty$ then we instead obtain
$\|\hat{\bm{b}}-\dot{\bm{b}}\|_{2}=O_{p}(m^{1/2}n_{U}^{-1/2})$. This implies
that, under the unconditional regime, a consistency result based on the
$l_{2}$-norm cannot hold for the entire vector of random effects when there is
a partnered fixed effect. If there is no partnered fixed effect though,
consistency of the entire vector is sometimes possible. For example, in the
Poisson counterexample, we demonstrate in the Appendix that
$\|\hat{\bm{b}}-\dot{\bm{b}}\|_{2}=O_{p}(m^{1/2}n^{-1/2})=o_{p}(1)$ when
$mn^{-1}\rightarrow 0$.

###  S2.5 Proof of Theorem 4 and Corollary 1

We will prove each of the three parts of the theorem separately. The proof of
part (a) also proves Corollary 1.

Part (a): When $mn_{U}^{-1}\rightarrow\infty$, we have from equation (S2.13)
and the order results for the remainder term in Section S3.2 that

$\displaystyle\bm{D}_{r}(\hat{\bm{b}}-\dot{\bm{b}})$
$\displaystyle=\bm{D}_{r}\bm{B}_{3}^{-1}\bm{D}_{r}\bm{D}_{r}^{-1}\dot{\phi}^{-1}\bm{Z}^{\top}(\bm{y}-\dot{\bm{\mu}})+o_{p}(1).$

This is identical to the proof of Theorem 3. Next, without loss of generality,
suppose $\bm{A}_{r}$ selects the first cluster only. Then we have

$\displaystyle{n}_{1}^{1/2}(\hat{\bm{b}}_{1}-\dot{\bm{b}}_{1})$
$\displaystyle=(n_{1}^{-1}\bm{X}_{1}^{\top}\dot{\bm{W}}_{1}\bm{X}_{1})^{-1}n_{1}^{-1/2}\\{\dot{\phi}^{-1}\bm{X}_{1}^{\top}(\bm{y}_{1}-\dot{\bm{\mu}}_{1})\\}+o_{p}(1)$
$\displaystyle\overset{\Delta}{=}\bm{P}_{n_{1}}+o_{p}(1).$

We wish to study the distribution of $\bm{P}_{n_{1}}$ as
$m,n_{L}\rightarrow\infty$. By definition,

$\displaystyle\underset{{m,n_{L}\rightarrow\infty}}{\lim}F_{\bm{P}_{n_{1}}}(\bm{x})=\underset{{m,n_{L}\rightarrow\infty}}{\lim}\int
F_{\bm{P}_{n_{1}}|\dot{\bm{b}}_{1}}(\bm{x})f(\dot{\bm{b}}_{1})d\dot{\bm{b}}_{1}.$

Since $F_{\bm{P}_{n_{1}}|\dot{\bm{b}}_{1}}(\bm{x})$ is a cdf, then
$F_{\bm{P}_{n_{1}}|\dot{\bm{b}}_{1}}(\bm{x})f(\dot{\bm{b}}_{1})$ is bounded by
$f(\dot{\bm{b}}_{1})$. Hence applying $\int
f(\dot{\bm{b}}_{1})d\dot{\bm{b}}_{1}=1$ and the dominated convergence theorem,
we obtain

$\displaystyle\underset{{m,n_{L}\rightarrow\infty}}{\lim}F_{\bm{P}_{n_{1}}}(\bm{x})$
$\displaystyle=\int\underset{{m,n_{L}\rightarrow\infty}}{\lim}F_{\bm{P}_{n_{1}}|\dot{\bm{b}}_{1}}(\bm{x})f(\dot{\bm{b}}_{1})d\dot{\bm{b}}_{1}=\int\Psi_{\bm{P}_{n_{1}}|\dot{\bm{b}}_{1}}(\bm{x})f(\dot{\bm{b}}_{1})d\dot{\bm{b}}_{1},$

where $\Psi_{\bm{P}_{n_{1}}|\dot{\bm{b}}_{1}}(\cdot)$ is the cdf associated
with $N(\bm{0},\bm{K}_{1})$, a result which follows from conditional
independence and the Lindeberg-Feller Central Limit Theorem used in Theorem 1.
The general result follows by noting that the same argument can be applied to
any finite subset of the random effects. Note also that the result holds
regardless of the true distribution of $\dot{\bm{b}}_{i}$.

Part (b): When $mn_{i}^{-1}\rightarrow\gamma_{i}\in(0,\infty)$, we have from
(S2.13) and the order results for the remainder term in Section S3.2 that

$\displaystyle{n}_{i}^{1/2}(\hat{\bm{b}}_{i}-\dot{\bm{b}}_{i})$
$\displaystyle=(n_{i}^{-1}\bm{X}_{i}^{\top}\dot{\bm{W}}_{i}\bm{X}_{i})^{-1}n_{i}^{-1/2}\\{\dot{\phi}^{-1}\bm{X}_{i}^{\top}(\bm{y}_{i}-\dot{\bm{\mu}}_{i})\\}-(\gamma_{i}m)^{-1/2}\sum_{i=1}^{m}\dot{\bm{b}}_{i}+O_{p}(n_{L}^{-1/2}),$

from the same development as in the proof of Part (a). Letting
$\bm{E}_{1}=(n_{i}^{-1}\bm{X}_{i}^{\top}\dot{\bm{W}}_{i}\bm{X}_{i})^{-1}n_{i}^{-1/2}\dot{\phi}^{-1}\bm{X}_{i}^{\top}(\bm{y}_{i}-\dot{\bm{\mu}}_{i})$
and $\bm{E}_{2}=m^{-1/2}\sum_{i=1}^{m}\dot{\bm{b}}_{i}$, then since
$\bm{E}_{1}$ and $\bm{E}_{2}$ are independent given $\dot{\bm{b}}_{i}$, we
obtain for any $i$,

$\displaystyle\lim_{m,n_{L}\rightarrow\infty}F_{\bm{E}_{1},\bm{E}_{2}}(\bm{x},\bm{y})$
$\displaystyle=\lim_{m,n_{L}\rightarrow\infty}\int
F_{\bm{E}_{1},\bm{E}_{2}|\dot{\bm{b}}_{i}}(\bm{x},\bm{y})f(\dot{\bm{b}}_{i})d\dot{\bm{b}}_{i}$
$\displaystyle=\lim_{m,n_{L}\rightarrow\infty}\int
F_{\bm{E}_{1}|\dot{\bm{b}}_{i}}(\bm{x})F_{\bm{E}_{2}|\dot{\bm{b}}_{i}}(\bm{y})f(\dot{\bm{b}}_{i})d\dot{\bm{b}}_{i}$
$\displaystyle=\int\lim_{m,n_{L}\rightarrow\infty}F_{\bm{E}_{1}|\dot{\bm{b}}_{i}}(\bm{x})F_{\bm{E}_{2}|\dot{\bm{b}}_{i}}(\bm{y})f(\dot{\bm{b}}_{i})d\dot{\bm{b}}_{i}$
$\displaystyle=\int\lim_{n_{L}\rightarrow\infty}F_{\bm{E}_{1}|\dot{\bm{b}}_{i}}(\bm{x})\lim_{m\rightarrow\infty}F_{\bm{E}_{2}|\dot{\bm{b}}_{i}}(\bm{y})f(\dot{\bm{b}}_{i})d\dot{\bm{b}}_{i}$
$\displaystyle=\Psi_{\bm{E}_{2}}(\bm{y})\int\lim_{n_{L}\rightarrow\infty}F_{\bm{E}_{1}|\dot{\bm{b}}_{i}}(\bm{x})f(\dot{\bm{b}}_{i})d\dot{\bm{b}}_{i},$

where $\Psi_{\bm{E}_{2}}(\cdot)$ is the cdf of $N(\bm{0},\dot{\bm{G}})$. The
third line follows from the Dominated Convergence Theorem since
$F_{\bm{E}_{1}|\dot{\bm{b}}_{i}}(\bm{x})$ and
$F_{\bm{E}_{2}|\dot{\bm{b}}_{i}}(\bm{y})$ are cdfs and $\int
f(\dot{\bm{b}}_{i})d\dot{\bm{b}}_{i}=1$. Thus $\bm{E}_{1}$ and $\bm{E}_{2}$
are asymptotically independent. The result follows from this asymptotic
independence.

Part (c): When $mn_{L}^{-1}\rightarrow 0$, we have from (S2.13) and the order
results for the remainder term in Section S3.2 that

$\displaystyle{m}^{1/2}(\hat{\bm{b}}-\dot{\bm{b}})$
$\displaystyle=-\bm{1}_{m}\otimes\bm{I}_{p}m^{-1/2}\sum_{i=1}^{m}\dot{\bm{b}}_{i}+o_{p}(1).$

The result then follows immediately from the normality assumption on
$\dot{\bm{b}}_{i}$.

###  S2.6 Proof of Theorem 5

Given $mn_{L}^{-2}\rightarrow 0$ and $mn_{U}^{-1/2}\rightarrow\infty$, by
summing equations (S2.12) and (S2.13) we see that the