# Tensor Generalized Estimating Equations
for Longitudinal Imaging Analysis

Xiang Zhang, Lexin Li, Hua Zhou, Dinggang Shen
and the Alzheimer’s Disease Neuroimaging Initiative

###### Abstract

In an increasing number of neuroimaging studies, brain images, which are in
the form of multidimensional arrays (tensors), have been collected on multiple
subjects at multiple time points. Of scientific interest is to analyze such
massive and complex longitudinal images to diagnose neurodegenerative
disorders and to identify disease relevant brain regions. In this article, we
treat those problems in a unifying regression framework with image predictors,
and propose tensor generalized estimating equations (GEE) for longitudinal
imaging analysis. The GEE approach takes into account intra-subject
correlation of responses, whereas a low rank tensor decomposition of the
coefficient array enables effective estimation and prediction with limited
sample size. We propose an efficient estimation algorithm, study the
asymptotics in both fixed $p$ and diverging $p$ regimes, and also investigate
tensor GEE with regularization that is particularly useful for region
selection. The efficacy of the proposed tensor GEE is demonstrated on both
simulated data and a real data set from the Alzheimer’s Disease Neuroimaging
Initiative (ADNI).

11footnotetext: Xiang Zhang is Graduate Student, Department of Statistics,
North Carolina State University, Raleigh, NC 27695-8203 (Email:
xzhang23@ncsu.edu). Lexin Li is Associate Professor, Division of
Biostatistics, University of California, Berkeley, Berkeley, CA 94720-3370
(Email: lexinli@berkeley.edu). Hua Zhou is Assistant Professor, Department of
Statistics, North Carolina State University, Raleigh, NC 27695-8203 (Email:
hua_zhou@ncsu.edu). Dinggang Shen is Professor, Department of Radiology,
University of North Carolina, Chapel Hill, NC 27599-7420 (E-mail:
dinggang_shen@med.unc.edu). The Alzheimer’s Disease Neuroimaging Initiative:
Data used in the preparation of this article were obtained from the ADNI data
base (http://adni.loni.usc.edu/). As such, the investigators within the ADNI
contributed to the design and implementation of ADNI and/or provided data but
did not participate in the analysis or writing of this article. A complete
listing of ADNI investigators is available at:
http://www.loni.usc.edu/ADNI/Data/ADNI_Authorship_List.pdf.

Key Words: Alzheimer’s disease; generalized estimating equations (GEE);
longitudinal imaging data; magnetic resonance imaging (MRI); multidimensional
array; tensor regression.

##  1 Introduction

Analyzing brain imaging data to study neuropsychiatric and neurodegenerative
disorders is gaining increasing interest in recent years (Lazar,, 2008;
Friston,, 2009; Hinrichs et al.,, 2009; Kang et al.,, 2012; Aston and Kirch,,
2012, among many others). There are a variety of forms, or modalities, of
images obtained through different imaging technologies, including magnetic
resonance imaging (MRI), functional magnetic resonance imaging (fMRI),
positron emission tomography (PET), and electroencephalography (EEG), among
others. Regardless of image modalities, it is of common scientific interest to
use brain images to diagnose neurodegenerative disorders, to predict onset of
neuropsychiatric diseases, and to identify disease relevant brain regions or
activity patterns. These problems can be collectively formulated as a
regression with a clinical outcome and an image predictor, whereas the image
takes a unifying form of multidimensional array, also known as _tensor_.

Early imaging studies typically involved only a handful of subjects. More
recently, a number of brain imaging databases are emerging with a relatively
large number of study subjects (ADHD,, 2014; ADNI,, 2014). Meanwhile, in an
increasing number of studies, images were acquired for each subject not only
at the baseline, but also over multiple visits, resulting in longitudinal
images. Our motivating example is a study from the Alzheimer’s Disease
Neuroimaging Initiative (ADNI). It consists of 88 subjects with mild cognitive
impairment (MCI), which is a prodromal stage of Alzheimer’s disease (AD). Each
subject had MRI scans at 5 different time points: baseline, 6-month, 12-month,
18-month and 24-month. After preprocessing, each MRI image is $32\times
32\times 32$ dimensional. Also measured for each subject at each visit was a
cognitive score, the Mini-Mental State Examination (MMSE), indicating
progression of the disease. It is scientifically important to understand
association between MCI/AD and the structural brain atrophy as reflected by
MRI. It is equally important to use MRI images to accurately predict AD/MCI,
as an accurate diagnosis is critical for timely therapy and possible delay of
the disease (Zhang et al.,, 2011).

While there has been an enormous literature on imaging analysis for AD, most
existing methods perform the prediction using only the baseline data, ignoring
data at the follow-up time points that often contain useful longitudinal
information. Recently, a small group of researchers started to use
longitudinal imaging data for individual-based classification (Misra et al.,,
2009; Davatzikos et al.,, 2009; McEvoy et al.,, 2011; Hinrichs et al.,, 2011)
and for cognitive score prediction (Zhang et al.,, 2012), whereas a limited
number of studies regressed longitudinal image responses on a collection of
covariates, first one voxel at a time then spatially smoothing the parameters
(Skup et al.,, 2012; Li et al.,, 2013). In general, longitudinal imaging
analysis is challenging, due to both the ultrahigh dimensionality and the
complex spatial structure of images, while the longitudinal correlation adds
another layer of complication.

Since the seminal work of Liang and Zeger, (1986), there has been a
substantive literature on statistical analysis of longitudinal data. See
Prentice and Zhao, (1991); Li, (1997); Qu et al., (2000); Xie and Yang,
(2003); Balan and Schiopu-Kratina, (2005); Song et al., (2009); Wang, (2011),
among many others. There is also a line of research studying variable
selection for longitudinal models, including Pan, (2001); Fan and Li, (2004);
Ni et al., (2010); Xue et al., (2010); Wang et al., (2012). However, all those
studies take the covariates as a vector, whereas in imaging regression,
covariates take the form of multi-dimensional arrays. Naively turning an array
into a vector would result in extremely high dimensionality. For instance, a
$32\times 32\times 32$ MRI image would require $32^{3}=32,768$ parameters.
Moreover, vectorization destroys inherent spatial information in images. There
have been some recent developments of statistical regression models for
image/tensor covariates; for instance, Caffo et al., (2010); Reiss and Ogden,
(2010); Wang et al., (2014). In particular, Zhou et al., (2013) proposed a
class of tensor regression models by imposing a low rank tensor decomposition
on the coefficient tensor. Although those methods directly work with a tensor
covariate, none has taken longitudinal tensors into account, and thus none is
immediately applicable to our longitudinal imaging study.

In this article, we propose tensor generalized estimating equations for
longitudinal imaging analysis. Our proposal consists of two key components: a
low rank tensor decomposition and generalized estimating equations (GEE).
Similar to Zhou et al., (2013), we choose to impose a low rank structure, the
CANDECOMP/PARAFAC (CP) decomposition (Kolda and Bader,, 2009), on the
coefficient array in GEE. This structure substantially reduces the number of
free parameters and makes subsequent estimation and inference feasible. But
unlike Zhou et al., (2013), we incorporate this low rank structure in
estimating equations to accommodate longitudinal correlation of the data. We
have chosen GEE over another popular approach, the mixed effects model, for
longitudinal imaging analysis. This is because the GEE approach only requires
the first two marginal moments and a working correlation structure for the
scalar response variable. By contrast, a mixed effects model requires
specification of a distribution for the parameters, which turns out to be a
formidable task for a tensor covariate. Within the tensor GEE framework, we
develop a scalable computation algorithm for solving the complicated tensor
estimating equations. Next we establish the asymptotic properties of the
solution of tensor GEE, including consistency and asymptotic normality under
two large sample scenarios: the number of parameters is fixed and the number
of parameters diverges along with the sample size. In particular, we show that
the tensor GEE estimator inherits the robustness feature of the classical GEE
estimator, in that the estimate is consistent even if the working correlation
structure is misspecified. Finally, we investigate regularization in the
context of tensor GEE. Regularization is crucial when the number of parameters
far exceeds the sample size, and is also useful for stabilizing estimates and
incorporating prior subject knowledge. For instance, employing an $L_{1}$
penalty in our tensor GEE in effect finds subregions of brains that are highly
relevant to the clinical outcome. This region selection is of scientific
interest itself, and corresponds to the intensively studied variable selection
problem in classical regressions with vector-valued predictors.

Our contributions are two-fold. First of all, our proposal offers a timely
response to the increasing availability of longitudinal imaging data along
with the growing interest of their analysis. To the best of our knowledge,
there has been very few systematic statistical methods developed for such an
analysis. Second, our work generalizes both the GEE approach from vector-
valued covariates to tensor-valued image covariate, as well as the tensor
regression model of Zhou et al., (2013) from independent imaging data to
longitudinal imaging data. Such a generalization parallels the extension in
classical regressions with vector predictors. This extension, however, is far
from trivial. Owing to the intrinsic complexity of both spatially and
temporally correlated observations as well as the huge data size, longitudinal
imaging analysis is much more challenging than both longitudinal analysis with
vector-valued predictors and imaging analysis at a single time point. Given
that the results of this kind are rare, our proposal offers a useful addition
to the literature of both longitudinal and imaging analysis.

The rest of the article is organized as follows. Section 2 proposes tensor GEE
for longitudinal imaging data, along with their estimation and regularization.
Section 3 presents the asymptotic results for the tensor GEE estimates.
Simulation studies and real data analysis are carried out in Sections 4 and 5,
respectively, followed by a discussion in Section 6.

##  2 Tensor Generalized Estimating Equations

###  2.1 Notations and Preliminaries

Suppose there are $n$ training subjects, and for the $i$-th subject, there are
observations over $m_{i}$ time points. For simplicity, we assume $m_{i}=m$ and
the time points are the same for all subjects. The observed data consist of
$\\{(Y_{ij},{\bm{X}}_{ij},{\bm{Z}}_{ij}),i=1,\ldots,n,j=1,\ldots,m\\}$, where,
for the $i$-th subject at the $j$-th time point, $Y_{ij}$ denotes the target
response, ${\bm{Z}}_{ij}\in\mathrm{I\\!R}\mathit{{}^{p_{0}}}$ is a
conventional predictor vector, and
${\bm{X}}_{ij}\in\mathrm{I\\!R}\mathit{{}^{p_{1}\times\dots\times p_{D}}}$ is
a $D$-dimensional array that represents the image covariate. The array
dimension $D$ depends on the image modality. With an image at a single time
point, for EEG, $D=2$, for MRI and PET, $D=3$, and for fMRI, $D=4$. Write
${\bm{Y}}_{i}=(Y_{i1},\dots,Y_{im})^{\mbox{\tiny{\sf T}}}$. A key attribute of
longitudinal data is that the observations from different subjects are
commonly assumed independent, but the observations from the same subject are
_correlated_. That is, the intra-subject covariance matrix,
$\textrm{Var}({\bm{Y}}_{i})\in\mathrm{I\\!R}\mathit{{}^{m\times m}}$, is not a
diagonal matrix but with some structure.

Next we review some key notations and operations of multidimensional array
that will be used throughout this article. The _inner product_ between two
tensors is defined as
$\langle{\bm{B}},{\bm{X}}\rangle=\langle\mathrm{vec}{\bm{B}},\mathrm{vec}{\bm{X}}\rangle=\sum_{i_{1},\ldots,i_{D}}\beta_{i_{1}\ldots
i_{D}}x_{i_{1}\ldots i_{D}}$, where the _$\mathrm{vec}({\bm{B}})$ operator_
stacks the entries of a tensor
${\bm{B}}\in\mathrm{I\\!R}\mathit{{}^{p_{1}\times\cdots\times p_{D}}}$ into a
column vector. The _outer product_ ,
${\bm{b}}_{1}\circ{\bm{b}}_{2}\circ\cdots\circ{\bm{b}}_{D}$, of $D$ vectors
${\bm{b}}_{d}\in\mathrm{I\\!R}\mathit{{}^{p_{d}}}$ is a
$p_{1}\times\cdots\times p_{D}$ array with entries
$({\bm{b}}_{1}\circ{\bm{b}}_{2}\circ\cdots\circ{\bm{b}}_{D})_{i_{1}\cdots
i_{D}}=\prod_{d=1}^{D}b_{di_{d}}$. The _mode- $d$ matricization_,
${\bm{B}}_{(d)}$, flattens a tensor ${\bm{B}}$ into a
$p_{d}\times\prod_{d^{\prime}\neq d}p_{d^{\prime}}$ matrix such that the
$(i_{1},\ldots,i_{D})$ element of the array ${\bm{B}}$ maps to the $(i_{d},j)$
element of the matrix ${\bm{B}}_{(d)}$, where $j=1+\sum_{d^{\prime}\neq
d}(i_{d^{\prime}}-1)\prod_{d^{\prime\prime}<d^{\prime},d^{\prime\prime}\neq
d}p_{d^{\prime\prime}}$.

A tensor ${\bm{B}}\in\mathrm{I\\!R}\mathit{{}^{p_{1}\times\cdots\times
p_{D}}}$ is said to admit a _rank- $R$ CANDECOMP/PARAFAC (CP) decomposition_
(Kolda and Bader,, 2009), if

$\displaystyle{\bm{B}}=\sum_{r=1}^{R}\mbox{\boldmath$\beta$}_{1}^{(r)}\circ\cdots\circ\mbox{\boldmath$\beta$}_{D}^{(r)},$
(1)

where
$\mbox{\boldmath$\beta$}_{d}^{(r)}\in\mathrm{I\\!R}\mathit{{}^{p_{d}}},d=1,\ldots,D,r=1,\ldots,R$,
are all column vectors, and ${\bm{B}}$ cannot be written as a sum of less than
$R$ outer products. The decomposition (1) is often represented by a shorthand,
${\bm{B}}=\llbracket{\bm{B}}_{1},\ldots,{\bm{B}}_{D}\rrbracket$, where
${\bm{B}}_{d}=[\mbox{\boldmath$\beta$}_{d}^{(1)},\ldots,\mbox{\boldmath$\beta$}_{d}^{(R)}]\in\mathrm{I\\!R}\mathit{{}^{p_{d}\times
R}}$. If a tensor
${\bm{B}}\in\mathrm{I\\!R}\mathit{{}^{p_{1}\times\cdots\times p_{D}}}$ admits
a rank-$R$ decomposition (1), then

$\displaystyle{\bm{B}}_{(d)}$
$\displaystyle={\bm{B}}_{d}({\bm{B}}_{D}\odot\cdots\odot{\bm{B}}_{d+1}\odot{\bm{B}}_{d-1}\odot\cdots\odot{\bm{B}}_{1})^{\mbox{\tiny{\sf
T}}}\,\text{
and}~{}~{}\mathrm{vec}\,{\bm{B}}=({\bm{B}}_{D}\odot\cdots\odot{\bm{B}}_{1}){\bf
1}_{R},$

where $\odot$ denotes the _Khatri-Rao product_ (Rao and Mitra,, 1971) of two
matrices ${\bm{B}}_{d}\in\mathrm{I\\!R}\mathit{{}^{p_{d}\times r}}$ and
${\bm{B}}_{d^{\prime}}\in\mathrm{I\\!R}\mathit{{}^{p_{d^{\prime}}\times r}}$
such that
${\bm{B}}_{d}\odot{\bm{B}}_{d^{\prime}}=\left[\mbox{\boldmath$\beta$}_{d}^{(1)}\otimes\mbox{\boldmath$\beta$}_{d^{\prime}}^{(1)}\;\mbox{\boldmath$\beta$}_{d}^{(2)}\otimes\mbox{\boldmath$\beta$}_{d^{\prime}}^{(2)}\;\ldots\;\mbox{\boldmath$\beta$}_{d}^{(R)}\otimes\mbox{\boldmath$\beta$}_{d^{\prime}}^{(R)}\right]\in\mathrm{I\\!R}\mathit{{}^{p_{d}p_{d^{\prime}}\times
R}}$, and $\otimes$ denotes the _Kronecker product_.

###  2.2 Tensor Generalized Estimating Equations

The GEE method has been widely employed for analyzing correlated longitudinal
data since the pioneer work of Liang and Zeger, (1986). It requires
specification of the first two moments of the conditional distribution of the
response given the covariates,
$\mu_{ij}=E(Y_{ij}|{\bm{X}}_{ij},{\bm{Z}}_{ij})$ and
$\sigma^{2}_{ij}=\text{Var}(Y_{ij}|{\bm{X}}_{ij},{\bm{Z}}_{ij})$. Following
Liang and Zeger, (1986), we assume $Y_{ij}$ is from an exponential family with
canonical link. Then

$\displaystyle\mu_{ij}({\bm{B}},\mbox{\boldmath$\gamma$})=\mu(\theta_{ij}),\;\textrm{
and
}\;\sigma^{2}_{ij}({\bm{B}},\mbox{\boldmath$\gamma$})=\phi\mu^{(1)}(\theta_{ij}),\quad
i=1,\ldots,n,\;j=1,\ldots,m,$

where $\mu(\cdot)$ is a differentiable canonical link function,
$\mu^{(1)}(\cdot)$ is its first derivative, $\theta_{ij}$ is the linear
systematic part, and $\phi$ is an over-dispersion parameter. In this article
we simply set $\phi=1$ while the extension to a general $\phi$ is
straightforward. $\theta_{ij}$ is associated with the covariates via the
relation

$\displaystyle\theta_{ij}=\mbox{\boldmath$\gamma$}^{\mbox{\tiny{\sf
T}}}{\bm{Z}}_{ij}+\langle{\bm{B}},{\bm{X}}_{ij}\rangle,$ (2)

where $\gamma$ is the coefficient vector associated with the covariate vector
${\bm{Z}}$, including the intercept, and ${\bm{B}}$ is the coefficient tensor
of the same size as ${\bm{X}}$ that captures effects of every array element of
${\bm{X}}$.

The GEE estimator of ${\bm{B}},\mbox{\boldmath$\gamma$}$ is then defined as
the solution of

$\displaystyle\sum_{i=1}^{n}\left\\{\frac{\partial\mbox{\boldmath$\mu$}_{i}({\bm{B}},\mbox{\boldmath$\gamma$})}{\partial[\mathrm{vec}({\bm{B}})^{\mbox{\tiny{\sf
T}}},\mbox{\boldmath$\gamma$}^{\mbox{\tiny{\sf T}}}]^{\mbox{\tiny{\sf
T}}}}\right\\}^{\mbox{\tiny{\sf
T}}}{\bm{V}}_{i}^{-1}\bigl{\\{}{\bm{Y}}_{i}-\mbox{\boldmath$\mu$}_{i}({\bm{B}},\mbox{\boldmath$\gamma$})\bigr{\\}}={\mathbf{0}},$
(3)

where ${\bm{Y}}_{i}=(Y_{i1},\dots,Y_{im})^{\mbox{\tiny{\sf T}}}$,
$\mbox{\boldmath$\mu$}_{i}({\bm{B}},\mbox{\boldmath$\gamma$})=[\mu_{i1}({\bm{B}},\mbox{\boldmath$\gamma$}),\ldots,\mu_{im}({\bm{B}},\mbox{\boldmath$\gamma$})]^{\mbox{\tiny{\sf
T}}}$, and ${\bm{V}}_{i}=\text{cov}({\bm{Y}}_{i})$ is the response covariance
matrix of the $i$-th subject. The first component in (3) is the derivative of
$\mbox{\boldmath$\mu$}_{i}({\bm{B}},\mbox{\boldmath$\gamma$})$ with respect to
the vector $[\mathrm{vec}({\bm{B}})^{\mbox{\tiny{\sf
T}}},\mbox{\boldmath$\gamma$}^{\mbox{\tiny{\sf T}}}]^{\mbox{\tiny{\sf
T}}}\in\mathrm{I\\!R}\mathit{{}^{p_{0}+\prod_{d}p_{d}}}$. As such, there are
totally $p_{0}+\prod_{d}p_{d}$ estimating equations to solve in (3). For
regression with image covariates, this dimension is ultrahigh and usually far
exceeds the sample size. For instance, for a regression with a $32\times
32\times 32$ MRI image predictor, an intercept, and two additional scalar
covariates, the number of equations to solve is in the scale of
$32^{3}+3=327,71$, resulting no unique solution when the sample size is only
in hundreds. It thus becomes crucial to reduce the number of estimating
equations.

Toward that end, we impose a low rank structure on the coefficient array
${\bm{B}}$. More specifically, we assume ${\bm{B}}$ in model (2) follows a CP
structure in (1),
${\bm{B}}=\llbracket{\bm{B}}_{1},\ldots,{\bm{B}}_{D}\rrbracket$, where
${\bm{B}}_{d}=[\mbox{\boldmath$\beta$}_{d}^{(1)},\ldots,\mbox{\boldmath$\beta$}_{d}^{(R)}]\in\mathrm{I\\!R}\mathit{{}^{p_{d}\times
R}}$. Then the systematic part in (2) becomes

$\displaystyle\theta_{ij}$ $\displaystyle=$
$\displaystyle\mbox{\boldmath$\gamma$}^{\mbox{\tiny{\sf
T}}}{\bm{Z}}_{ij}+\langle\sum_{r=1}^{R}\mbox{\boldmath$\beta$}_{1}^{(r)}\circ\cdots\circ\mbox{\boldmath$\beta$}_{D}^{(r)},{\bm{X}}_{ij}\rangle$
(4) $\displaystyle=$ $\displaystyle\mbox{\boldmath$\gamma$}^{\mbox{\tiny{\sf
T}}}{\bm{Z}}_{ij}+\langle({\bm{B}}_{D}\odot\cdots\odot{\bm{B}}_{1}){\bf
1}_{R},\mathrm{vec}{\bm{X}}_{ij}\rangle.$

Adopting (4), we propose the _tensor generalized estimating equations_
estimator of ${\bm{B}},\mbox{\boldmath$\gamma$}$, defined as the solution of

$\displaystyle\sum_{i=1}^{n}\left\\{\frac{\partial\mbox{\boldmath$\mu$}_{i}({\bm{B}},\mbox{\boldmath$\gamma$})}{\partial[\mbox{\boldmath$\beta$}_{\bm{B}}^{\mbox{\tiny{\sf
T}}},\mbox{\boldmath$\gamma$}^{\mbox{\tiny{\sf T}}}]^{\mbox{\tiny{\sf
T}}}}\right\\}^{\mbox{\tiny{\sf
T}}}{\bm{V}}_{i}^{-1}\bigl{\\{}{\bm{Y}}_{i}-\mbox{\boldmath$\mu$}_{i}({\bm{B}},\mbox{\boldmath$\gamma$})\bigr{\\}}={\mathbf{0}},$
(5)

where
$\mbox{\boldmath$\beta$}_{\bm{B}}=\mathrm{vec}({\bm{B}}_{1},\ldots,{\bm{B}}_{D})$,
and the subscript ${\bm{B}}$ is to remind that $\beta$ is constructed based on
the CP decomposition of a given coefficient tensor
${\bm{B}}=\llbracket{\bm{B}}_{1},\ldots,{\bm{B}}_{D}\rrbracket$. Comparing to
the classical GEE (3), the derivative is now with respect to
$[\mbox{\boldmath$\beta$}_{\bm{B}}^{\mbox{\tiny{\sf
T}}},\mbox{\boldmath$\gamma$}^{\mbox{\tiny{\sf T}}}]^{\mbox{\tiny{\sf
T}}}\in\mathrm{I\\!R}\mathit{{}^{p_{0}+R\sum_{d}p_{d}}}$. Consequently, the
number of estimating equations has reduced from the exponential order
$p_{0}+\prod_{d}p_{d}$ to the linear order $p_{0}+R\sum_{d}p_{d}$. This
substantial reduction in dimensionality, as we will demonstrate later, enables
effective estimation and inference, and also provides a sound recovery of both
low rank and high rank signals.

Examining (5), the true intra-subject covariance structure ${\bm{V}}_{i}$ is
usually unknown in practice. The classical GEE adopts a working covariance
matrix, specified through a working correlation matrix ${\bm{R}}$. That is,
${\bm{V}}_{i}={\bm{A}}^{1/2}_{i}({\bm{B}},\mbox{\boldmath$\gamma$}){\bm{R}}{\bm{A}}^{1/2}_{i}({\bm{B}},\mbox{\boldmath$\gamma$})$,
where ${\bm{A}}_{i}({\bm{B}},\mbox{\boldmath$\gamma$})$ is an $m\times m$
diagonal matrix with $\sigma^{2}_{ij}({\bm{B}},\mbox{\boldmath$\gamma$})$ on
the diagonal and ${\bm{R}}$ is the $m$-by-$m$ working intra-subject
correlation matrix. Some commonly used correlation structures include
independence, autocorrelation (AR), compound symmetry, and unstructured
correlation, among others. The correlation matrix ${\bm{R}}$ may involve
additional parameters, which can be estimated using residual-based moment
method.

By both adopting this working covariance/correlation idea, and explicitly
evaluating the derivative in (5), we finally arrive at the formal definition
of the tensor GEE estimator, which is the solution
$(\widehat{\bm{B}},\hat{}\mbox{\boldmath$\gamma$})$ of the following
estimating equations

$\sum_{i=1}^{n}\begin{pmatrix}[{\bm{J}}_{1}\ldots{\bm{J}}_{D}]^{\mbox{\tiny{\sf
T}}}\mathrm{vec}({\bm{X}}_{i})\\\
({\bm{Z}}_{i1},\ldots,{\bm{Z}}_{im})\end{pmatrix}{\bm{A}}^{1/2}_{i}({\bm{B}},\mbox{\boldmath$\gamma$}){\bf\widehat{R}}^{-1}{{\bm{A}}}^{-1/2}_{i}({\bm{B}},\mbox{\boldmath$\gamma$})\bigl{\\{}{\bm{Y}}_{i}-\mbox{\boldmath$\mu$}_{i}({\bm{B}},\mbox{\boldmath$\gamma$})\bigr{\\}}={\mathbf{0}},$
(6)

where $\widehat{{\bm{R}}}$ is an estimated correlation matrix,
$\mathrm{vec}({\bm{X}}_{i})=(\mathrm{vec}({\bm{X}}_{i1}),\ldots,\mathrm{vec}({\bm{X}}_{im}))$
is a $\prod_{d=1}^{D}p_{d}\times m$ matrix, ${\bm{J}}_{d}$ is the
$\prod_{d=1}^{D}p_{d}\times p_{d}R$ Jacobian matrix of the form
$\mbox{\boldmath$\Pi$}_{d}[({\bm{B}}_{D}\odot\cdots\odot{\bm{B}}_{d+1}\odot{\bm{B}}_{d-1}\odot\cdots\odot{\bm{B}}_{1})\otimes{\bm{I}}_{p_{d}}]$,
where $\mbox{\boldmath$\Pi$}_{d}$ is the
$(\prod_{d=1}^{D}p_{d})$-by-$(\prod_{d=1}^{D}p_{d})$ permutation matrix that
reorders $\mathrm{vec}{\bm{B}}_{(d)}$ to obtain $\mathrm{vec}{\bm{B}}$, i.e.,
$\mathrm{vec}{\bm{B}}=\mbox{\boldmath$\Pi$}_{d}\,\mathrm{vec}{\bm{B}}_{(d)}$.
Note that $\mu^{(1)}(\theta_{ij})$ has been canceled by the diagonal of the
matrix ${\bm{A}}^{-1}_{i}$ due to the property of canonical link. For ease of
presentation, we denote the left hand side of equation (6) as
${\bm{s}}({\bm{B}},\mbox{\boldmath$\gamma$})$, and write the tensor GEE (6) as
${\bm{s}}({\bm{B}},\mbox{\boldmath$\gamma$})={\mathbf{0}}$.

###  2.3 Estimation

Directly solving the tensor generalized estimating equations (6) with respect
to $({\bm{B}},\mbox{\boldmath$\gamma$})$ can be computational intensive, as
the mean function of the response given the covariates is nonlinear in the
parameters and the Jacobian matrices ${\bm{J}}_{1},\ldots,{\bm{J}}_{D}$ also
depend on the unknown parameters. We propose to iteratively solve the sub-GEE
for ${\bm{B}}_{1},\ldots,{\bm{B}}_{D}$, along with $\gamma$, one at a time,
while keeping all other components fixed. When updating
${\bm{B}}_{d}\in\mathrm{I\\!R}\mathit{{}^{p_{d}\times R}}$, the systematic
part $\theta_{ij}({\bm{B}},\mbox{\boldmath$\gamma$})$ can be rewritten as

$\displaystyle\theta_{ij}({\bm{B}},\mbox{\boldmath$\gamma$})$ $\displaystyle=$
$\displaystyle\mbox{\boldmath$\gamma$}^{\mbox{\tiny{\sf
T}}}{\bm{Z}}_{ij}+\langle{\bm{B}},{\bm{X}}_{ij}\rangle$ $\displaystyle=$
$\displaystyle\mbox{\boldmath$\gamma$}^{\mbox{\tiny{\sf
T}}}{\bm{Z}}_{ij}+\langle{\bm{B}}_{d},{\bm{X}}_{ij(d)}({\bm{B}}_{D}\odot\cdots\odot{\bm{B}}_{d+1}\odot{\bm{B}}_{d-1}\odot\cdots\odot{\bm{B}}_{1})\rangle,$

where ${\bm{X}}_{ij(d)}$ is the mode-$d$ matricization of the tensor
${\bm{X}}_{ij}$. As such, the systematic part
$\theta_{ij}({\bm{B}},\mbox{\boldmath$\gamma$})$ becomes linear in
${\bm{B}}_{d}$. The Jacobian matrix ${\bm{J}}_{d}$ is free of ${\bm{B}}_{d}$
and depends on the covariates and fixed parameters only. Consequently, each
step reduces to a standard GEE problem with $Rp_{d}$ parameters, which can be
solved using standard statistical softwares.

A problem of practical interest is to choose the rank $R$ for ${\bm{B}}$ in
its CP decomposition. This can be viewed as a _model selection_ problem. Pan,
(2001) proposed a quasi-likelihood independence model criterion for the
classical GEE model selection, by evaluating the likelihood under the
independence working correlation assumption. In our tensor GEE setup, we use
the following BIC-type information criterion

$\text{BIC}(R)=-2\ell(\widehat{\bm{B}}(R),\hat{}\mbox{\boldmath$\gamma$};{\bm{I}}_{m})+\log(n)p_{e},$
(7)

where $\ell(\widehat{\bm{B}}(R),\hat{}\mbox{\boldmath$\gamma$};{\bm{I}}_{m})$
is the log-likelihood evaluated at the tensor GEE estimator
$\hat{}\mbox{\boldmath$\gamma$}$ and $\widehat{\bm{B}}(R)$ with a working rank
$R$ and the independence working correlation structure ${\bm{I}}_{m}$. For
simplicity, we call this criterion BIC, as the term $\log(n)$ is used. Because
the CP decomposition itself is not unique, but can be made so under some minor
conditions (Zhou et al.,, 2013), the actual number of estimating equations, or
the effective number of parameters, is of the form:
$p_{e}=R(p_{1}+p_{2})-R^{2}$ for $D=2$, and $p_{e}=R(\sum_{d}p_{d}-D+1)$ for
$D>2$. We choose $R$ that minimizes this criterion among a series of working
ranks. We will briefly illustrate its use in Section 4.1.

###  2.4 Regularization

Even after introducing a low rank structure in our tensor GEE, regularization
can still be useful, as the number of subjects is often limited in a
neuroimaging study. In this section, we consider a general form of regularized
tensor GEE that includes a variety of penalty functions. Then in Section 4.3,
we will illustrate with a lasso penalty that is capable of identifying sub-
regions of brains associated with the clinical outcome. Specifically, we
consider the following regularized tensor GEE

$\displaystyle{\bm{s}}({\bm{B}},\mbox{\boldmath$\gamma$})+\begin{pmatrix}\partial_{\beta_{11}^{(1)}}P_{\lambda}(|\beta_{11}^{(1)}|,\rho)\\\
\vdots\\\ \partial_{\beta_{di}^{(r)}}P_{\lambda}(|\beta_{di}^{(r)}|,\rho)\\\
\vdots\\\
\partial_{\beta_{Dp_{D}}^{(R)}}P_{\lambda}(|\beta_{Dp_{D}}^{(R)}|,\rho)\end{pmatrix}={\mathbf{0}}_{p_{e}},$

where $P_{\lambda}(|\beta|,\rho)$ is a scalar penalty function, $\rho$ is the
penalty tuning parameter, $\lambda$ is an index for the penalty family,
$\partial_{\beta}P_{\lambda}(|\beta|,\rho)$ is the subgradient with respect to
argument $\beta$, and the subscript $p_{e}$ of ${\mathbf{0}}$ is a reminder of
the number of estimating equations to solve. Some widely used penalties
include: power family (Frank and Friedman,, 1993), in which
$P_{\lambda}(|\beta|,\rho)=\rho|\beta|^{\lambda}$, $\lambda\in(0,2]$, and in
particular lasso (Tibshirani,, 1996) ($\lambda=1$) and ridge ($\lambda=2$);
elastic net (Zou and Hastie,, 2005), in which
$P_{\lambda}(|\beta|,\rho)=\rho[(\lambda-1)\beta^{2}/2+(2-\lambda)|\beta|],\lambda\in[1,2]$;
and SCAD (Fan and Li,, 2001), in which
$\partial/\partial|\beta|P_{\lambda}(|\beta|,\rho)=\rho\left\\{1_{\\{|\beta|\leq\rho\\}}+(\lambda\rho-|\beta|)_{+}/(\lambda-1)\rho
1_{\\{|\beta|>\rho\\}}\right\\}$, $\lambda>2$, among many others.

Thanks to the separability of parameters in the regularization term, the
alternating updating strategy still applies. When updating ${\bm{B}}_{d}$, we
solve the penalized sub-GEE

$\displaystyle{\bm{s}}_{d}({\bm{B}}_{d})+\begin{pmatrix}\partial_{\mbox{\boldmath$\beta$}_{d1}^{(1)}}P_{\lambda}(|\beta_{d1}^{(1)}|,\rho)\\\
\vdots\\\
\partial_{\mbox{\boldmath$\beta$}_{di}^{(r)}}P_{\lambda}(|\beta_{di}^{(r)}|,\rho)\\\
\vdots\\\
\partial_{\mbox{\boldmath$\beta$}_{dp_{d}}^{(R)}}P_{\lambda}(|\beta_{dp_{D}}^{(R)}|,\rho)\end{pmatrix}={\mathbf{0}}_{Rp_{d}},$
(8)

where ${\bm{s}}_{d}$ is the sub-estimation equation for block ${\bm{B}}_{d}$,
and there are $Rp_{d}$ equations to solve at this step. Anti-derivative of
${\bm{s}}_{d}$ is recognized as the loss of an Aitken linear model with block
diagonal covariance matrix. Thus after linear transformation of ${\bm{Y}}_{i}$
and the working design matrix, solution to (8) is same as the minimizer of a
regular penalized weighted least squares problem, for which many software
packages exist. The fitting procedure boils down to alternating penalized
weighted least squares problem.

##  3 Theory

In this section, we study the asymptotic properties of the unregularized
tensor GEE estimator as the number of subjects $n$ goes to infinity, while we
assume the true rank of the tensor coefficient is known. We investigate two
scenarios: the number of parameters is fixed in Section 3.1, and the number of
parameters diverges in Section 3.2. For ease of exposition, we omit the
vector-valued covariates ${\bm{Z}}$ and the associated parameters $\gamma$,
while the results can be easily extended to incorporate them. Our development
builds upon and extends the previous work of Xie and Yang, (2003); Balan and
Schiopu-Kratina, (2005); Wang, (2011) from classical vector GEE to tensor GEE,
while we spell out the similarity as well as difference in asymptotics when
comparing the vector and tensor GEE. We show that tensor GEE estimator
inherits the key advantage of the classical GEE estimator in that it remains
consistent even if the working correlation structure is misspecified. On the
other hand, we note that, although one can generalize the classical GEE
asymptotics by directly vectorizing the tensor, it would have to require a
more stringent set of conditions. By contrast, we could achieve the robustness
in consistency for our tensor GEE based on a weaker set of conditions, and we
achieve this by imposing and exploiting the special structure of the
coefficient tensor.

###  3.1 Asymptotics for Fixed Dimension

We begin with the list of regularity conditions for the asymptotics of tensor
GEE with a fixed number of parameters.

  1. (A1)

The elements of ${\bm{X}}_{ij}$, $i=1,\dots,n$, $j=1,\dots,m$, are uniformly
bounded by a finite constant.

  2. (A2)

The true value ${\bm{B}}_{0}$ of the unknown parameter lies in the interior of
a compact parameter space $\mathcal{B}$ and follows a rank-$R$ CP structure
defined in (1).

  3. (A3)

Letting
$I({\bm{B}})=n^{-1}\sum_{i=1}^{n}[{\bm{J}}_{1}\ldots{\bm{J}}_{D}]^{\mbox{\tiny{\sf
T}}}\mathrm{vec}({\bm{X}}_{i})\mathrm{vec}({\bm{X}}_{i})^{\mbox{\tiny{\sf
T}}}[{\bm{J}}_{1}\ldots{\bm{J}}_{D}]$. It is assumed that there exist two
positive constants $c_{1}<c_{2}$ such that

$c_{1}\leq\lambda_{\min}(I({\bm{B}}))\leq\lambda_{\max}(I({\bm{B}}))\leq
c_{2},$

over the set
$\\{{\bm{B}}:||\mbox{\boldmath$\beta$}_{{\bm{B}}}-\mbox{\boldmath$\beta$}_{{\bm{B}}_{0}}||\leq\triangle
n^{-1/2}\\}$ for some constant $\triangle>0$, where $\lambda_{\min}$ and
$\lambda_{\max}$ are smallest and largest eigenvalue, respectively. It is also
assumed that on the same set $I({\bm{B}})$ has a constant rank.

  4. (A4)

The true intra-subject correlation matrix ${\bm{R}}_{0}$ has bounded
eigenvalues from zero and infinity. The estimated working correlation matrix
satisfies
$\|\widehat{{\bm{R}}}^{-1}-\tilde{{\bm{R}}}^{-1}\|_{\textbf{F}}=O_{p}(n^{-1/2})$,
where $\|\cdot\|_{\textbf{F}}$ is the Frobenius norm, $\tilde{{\bm{R}}}$ is
some positive definite matrix with bounded eigenvalues from zero and infinity,
and $\tilde{{\bm{R}}}={\bm{R}}_{0}$ is _not_ required.

  5. (A5)

For some constant $\delta>0$ and $M_{1}>0$,
$E(\|{\bm{A}}^{-1/2}_{i}({\bm{B}}_{0})({\bm{Y}}_{i}-\mbox{\boldmath$\mu$}_{i}({\bm{B}}_{0}))\|)^{2+\delta}\leq
M_{1}$ for all $1\leq i\leq n$, where ${\bm{A}}^{-1/2}_{i}({\bm{B}}_{0})$ is
the covariance matrix of ${\bm{Y}}_{i}$.

  6. (A6)

$\sigma^{-1}_{ij}({\bm{B}}_{0})(Y_{ij}-\mu_{ij}({\bm{B}}_{0}))$ has sub-
Gaussian tails for all $i=1,\dots,n$, $j=1,\dots,m$.

  7. (A7)

The elements of
$\partial\theta_{ij}(\mbox{\boldmath$\beta$}_{{\bm{B}}_{0}})/\partial\mbox{\boldmath$\beta$}_{{\bm{B}}_{0}}$,
$i=1,\dots,n$, $j=1,\dots,m$, are uniformly bounded by a finite constant.

  8. (A8)

Denote $\mu^{(k)}(\theta_{ij})$ the $k$-th derivative of $\mu(\theta_{ij})$,
where $\theta_{ij}$ is the linear systematic part evaluated at the GEE
solution $\widehat{\bm{B}}$. It is assumed that $\mu^{(1)}(\theta_{ij})$ are
uniformly bounded away from zero and infinity, and $\mu^{(k)}(\theta_{ij})$
are uniformly bounded by a finite constant, over the set
$\\{{\bm{B}}:||\mbox{\boldmath$\beta$}_{{\bm{B}}}-\mbox{\boldmath$\beta$}_{{\bm{B}}_{0}}||\leq\triangle
n^{-1/2}\\}$, for some constant $\triangle>0$, $i=1,\dots,n$, $j=1,\dots,m$,
and $k=2,3$.

  9. (A9)

Denote
${\bm{H}}({\bm{B}},{\bm{X}}_{ij})=\partial[{\bm{J}}_{1}\cdots{\bm{J}}_{D}]^{\mbox{\tiny{\sf
T}}}\text{vec}({\bm{X}}_{ij})/\partial\text{vec}^{\mbox{\tiny{\sf
T}}}({\bm{B}})$. ${\bm{H}}({\bm{B}},{\bm{X}}_{ij})$ is the Hessian of the
linear systematic part $\theta_{ij}$ under tensor structure. There exist two
positive constants $c_{3}<c_{4}$ such that

$c_{3}\leq\lambda_{\min}({\bm{H}}({\bm{B}},{\bm{X}}_{ij}))\leq\lambda_{\max}({\bm{H}}({\bm{B}},{\bm{X}}_{ij}))\leq
c_{4},$

over the set
$\\{{\bm{B}}:||\mbox{\boldmath$\beta$}_{{\bm{B}}}-\mbox{\boldmath$\beta$}_{{\bm{B}}_{0}}||\leq\triangle
n^{-1/2}\\}$ for some constant $\triangle>0$, $i=1,\dots,n$ and $j=1,\dots,m$.

A few remarks are in order. Conditions (A2) and (A3) are required for model
identifiability of tensor GEE (Zhou et al.,, 2013). We observe that, the
matrix $I({\bm{B}})$ in (A3) is an $R\sum_{d=1}^{D}p_{d}\times
R\sum_{d=1}^{D}p_{d}$ matrix, and thus (A3) is much weaker than the
nonsingularity condition on the design matrix if one were to directly
vectorize the tensor covariate. Condition (A4) is commonly imposed in the GEE
literature. It only requires a consistent estimator $\widehat{{\bm{R}}}$ of
some $\tilde{\bm{R}}$, in the sense
$\|\widehat{{\bm{R}}}^{-1}-\tilde{{\bm{R}}}^{-1}\|_{\textbf{F}}=O_{p}(n^{-1/2})$.
$\tilde{\bm{R}}$ needs to be well behaved in that it is positive definite with
bounded eigenvalues from zero and infinity, but $\tilde{\bm{R}}$ does _not_
have to be the true intra-subject correlation ${\bm{R}}$. This condition
essentially leads to the robust feature in Theorem 1 that the tensor GEE
estimate is consistent even if the working correlation structure is
misspecified. Conditions (A5) and (A6) regulate the tail behavior of the
residuals so that the noise cannot accumulate too fast, and we can employ the
Lindeberg-Feller central limit theorem to control the asymptotic behavior of
the residuals. Condition (A7) states the gradients of the systematic part
evaluated at the truth are well-defined. Condition (A8) concerns the canonical
link and generally holds for common exponential families, for example, the
binomial distribution with
$\mu(\theta_{ij})=\exp{\theta_{ij}}/(1+\exp{\theta_{ij}})$, and the Poisson
distribution with $\mu(\theta_{ij})=\exp{\theta_{ij}}$. Condition (A9) ensures
that the Hessian matrix of the linear systematic part, which is highly sparse,
is well-behaved in a neighborhood of the true value.

Before we turn to the asymptotics of the tensor GEE estimator, we address two
components involved in the estimating equations: the initial estimator and the
correlation estimator. Recall the tensor GEE estimator $\widehat{\bm{B}}$ is
obtained by solving the equations in (6). After dropping the covariate vector
${\bm{Z}}$, the tensor estimating equations become

$\sum_{i=1}^{n}[{\bm{J}}_{1}\ldots{\bm{J}}_{D}]^{\mbox{\tiny{\sf
T}}}\mathrm{vec}({\bm{X}}_{i}){\bm{A}}^{1/2}_{i}({\bm{B}})\widehat{{\bm{R}}}^{-1}{\bm{A}}^{-1/2}_{i}({\bm{B}})\bigl{\\{}{\bm{Y}}_{i}-\mbox{\boldmath$\mu$}_{i}({\bm{B}})\bigr{\\}}={\mathbf{0}},$
(9)

where $\widehat{{\bm{R}}}$ is any estimator of the intra-subject correlation
matrix satisfying the condition (A4). We still denote the left hand side by
${\bm{s}}({\bm{B}})$. Note that (9) involves the unknown correlation
${\bm{R}}$, and its estimate $\widehat{{\bm{R}}}$ is often obtained via
residual-based moment method, which in turn requires an initial estimator of
${\bm{B}}$. Next, we examine some frequently used estimators of
$\widehat{{\bm{B}}}$ and $\widehat{{\bm{R}}}$.

A customary initial estimator $\widehat{{\bm{B}}}$ in the GEE literature is
the one that assumes an independent working correlation. That is, one
completely ignores possible intra-subject correlation, and the corresponding
tensor GEE becomes

$\sum_{i=1}^{n}[{\bm{J}}_{1}\ldots{\bm{J}}_{D}]^{\mbox{\tiny{\sf
T}}}\mathrm{vec}({\bm{X}}_{i})\bigl{\\{}{\bm{Y}}_{i}-\mbox{\boldmath$\mu$}_{i}({\bm{B}})\bigr{\\}}={\mathbf{0}}.$

Denoting the equations as ${\bm{s}}_{init}({\bm{B}})={\mathbf{0}}$, and the
solution as $\widehat{\bm{B}}_{init}$, the next Lemma shows that it is a
consistent estimator of the true ${\bm{B}}_{0}$.

######  Lemma 1.

Under conditions (A1)-(A3) and (A5)-(A9), there exists a root
$\widehat{\bm{B}}_{init}$ of the equations
${\bm{s}}_{init}({\bm{B}})={\mathbf{0}}$ satisfing that

$\displaystyle\|\mbox{\boldmath$\beta$}_{\widehat{{\bm{B}}}_{init}}-\mbox{\boldmath$\beta$}_{{\bm{B}}_{0}}\|=O_{p}(n^{-1/2}).$

Here
$\mbox{\boldmath$\beta$}_{\bm{B}}=\mathrm{vec}({\bm{B}}_{1},\ldots,{\bm{B}}_{D})$,
and is constructed based on the CP decomposition of a given tensor
${\bm{B}}=\llbracket{\bm{B}}_{1},\ldots,{\bm{B}}_{D}\rrbracket$, as defined
before.

Given a consistent initial estimator of ${\bm{B}}_{0}$, there exist multiple
choices for the working correlation structure, e.g., autocorrelation, compound
symmetry, and the nonparametric structure (Balan and Schiopu-Kratina,, 2005).
We will investigate those choices in our simulations and real data analysis.

Next we establish the consistency and asymptotic normality of the tensor GEE
estimator from (9).

######  Theorem 1.

Under conditions (A1)-(A9), there exists a root $\widehat{\bm{B}}$ of the
equations ${\bm{s}}({\bm{B}})={\mathbf{0}}$ satisfing that

$\displaystyle\|\mbox{\boldmath$\beta$}_{\widehat{{\bm{B}}}}-\mbox{\boldmath$\beta$}_{{\bm{B}}_{0}}\|=O_{p}(n^{-1/2}).$

The key message of Theorem 1, as implied by condition (A4), is that the
consistency of the tensor coefficient estimator $\widehat{\bm{B}}$ does _not_
require the estimated working correlation $\widehat{{\bm{R}}}$ being a
consistent estimator of the true correlation ${\bm{R}}$. This protects us from
potential misspecification of the intra-subject correlation structure. Such a
robustness feature is well known for GEE estimator with vector-valued
covariates. Theorem 1 confirms and extends this result to the tensor GEE case
with image covariates. We also remark that, although the asymptotics of the
classical GEE can in principle be generalized to the tensor data by directly
vectorizing the coefficient array, the ultrahigh dimensionality of the
parameters would have made the regularity conditions such as (A3) unrealistic.
By contrast, Theorem 1 ensures that one could still enjoy the consistency and
robustness properties, by taking into account the structural information of
the tensor coefficient under the GEE framework.

Under condition (A4), we define

$\displaystyle\tilde{\bm{M}}_{n}({\bm{B}})$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{n}[{\bm{J}}_{1}\ldots{\bm{J}}_{D}]^{\mbox{\tiny{\sf
T}}}\mathrm{vec}({\bm{X}}_{i}){\bm{A}}^{1/2}_{i}({\bm{B}})\tilde{{\bm{R}}}^{-1}{\bm{R}}_{0}\tilde{{\bm{R}}}^{-1}{\bm{A}}^{1/2}_{i}({\bm{B}})\mathrm{vec}^{\mbox{\tiny{\sf
T}}}({\bm{X}}_{i})[{\bm{J}}_{1}\ldots{\bm{J}}_{D}],$
$\displaystyle\tilde{\bm{D}}_{n1}({\bm{B}})$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{n}[{\bm{J}}_{1}\ldots{\bm{J}}_{D}]^{\mbox{\tiny{\sf
T}}}\mathrm{vec}({\bm{X}}_{i}){\bm{A}}^{1/2}_{i}({\bm{B}})\tilde{{\bm{R}}}^{-1}{\bm{A}}^{1/2}_{i}({\bm{B}})\mathrm{vec}^{\mbox{\tiny{\sf
T}}}({\bm{X}}_{i})[{\bm{J}}_{1}\ldots{\bm{J}}_{D}].$

As we will show in the appendix, $\tilde{\bm{M}}_{n}({\bm{B}})$ approximates
the covariance matrix of ${\bm{s}}({\bm{B}})$ in (9), while
$\tilde{\bm{D}}_{n1}({\bm{B}})$ approximates the leading term of the negative
gradient of ${\bm{s}}({\bm{B}})$ with respect to
$\mbox{\boldmath$\beta$}_{\bm{B}}$. Then the next theorem gives the asymptotic
normality of the tensor GEE estimator.

######  Theorem 2.

Under conditions (A1)-(A9), for any vector
${\bm{b}}\in\mathrm{I\\!R}\mathit{{}^{R\sum_{d=1}^{D}p_{d}}}$ such that
$\|{\bm{b}}\|=1$, we have

$\displaystyle{\bm{b}}^{\mbox{\tiny{\sf
T}}}\tilde{{\bm{M}}_{n}}^{-1/2}({\bm{B}}_{0})\tilde{{\bm{D}}}_{n1}({\bm{B}}_{0})\bigl{(}\mbox{\boldmath$\beta$}_{\widehat{{\bm{B}}}}-\mbox{\boldmath$\beta$}_{{\bm{B}}_{0}}\bigr{)}\to\mathrm{Normal}(0,1)\textrm{
in distribution.}$

By Theorem 2 and Cramér-Wold theorem, one can derive the sandwich covariance
estimator of $\textrm{Var}(\mbox{\boldmath$\beta$}_{\widehat{\bm{B}}})$, and
carry out the subsequent Wald inference. Specifically, it is easy to see that
the variance of the GEE estimator can be approximated by the asymptotic
variance
$\tilde{\bm{D}}_{n1}^{-1}({\bm{B}}_{0})\tilde{\bm{M}}_{n}({\bm{B}}_{0})\tilde{\bm{D}}_{n1}^{-1}({\bm{B}}_{0})$.
Since it involves the unknown terms ${\bm{B}}_{0},{\bm{R}}_{0}$ and
$\tilde{\bm{R}}$, we plug in, respectively, $\widehat{{\bm{B}}}$,
$n^{-1}\sum_{i=1}^{n}{\bm{A}}^{-1/2}_{i}(\widehat{{\bm{B}}})\\{{\bm{Y}}_{i}-\mbox{\boldmath$\mu$}_{i}(\widehat{{\bm{B}}})\\}\\{{\bm{Y}}_{i}-\mbox{\boldmath$\mu$}_{i}(\widehat{{\bm{B}}})\\}^{\mbox{\tiny{\sf
T}}}{\bm{A}}^{-1/2}_{i}(\widehat{{\bm{B}}})$, and $\widehat{{\bm{R}}}$, which
leads to the sandwich estimator,

$\displaystyle\widehat{\text{Var}}(\mbox{\boldmath$\beta$}_{\widehat{{\bm{B}}}})=\tilde{\bm{D}}_{n1}^{-1}(\widehat{{\bm{B}}})\tilde{\bm{M}}_{n}(\widehat{{\bm{B}}})\tilde{\bm{D}}_{n1}^{-1}(\widehat{{\bm{B}}}).$

This sandwich formula in turn can be used to construct asymptotic confidence
interval or asymptotic hypothesis testing through the usual Wald inference.

###  3.2 Asymptotics for Diverging Dimension

We next study the asymptotics when the number of parameters diverges. We
assume that $p_{d}\sim p_{n}$ for $d=1,\dots,D$, where $a_{n}\sim b_{n}$ means
$a_{n}=O(b_{n})$ and $b_{n}=O(a_{n})$. We also assume that the rank $R$ is
fixed in the tensor GEE. Next we list the required regularity conditions.
Since the conditions (A1), (A2), (A5)–(A7) are the same as in Section 3.1, we
only list the conditions that are different, while we relabel those same
conditions as (A1∗), (A2∗), (A5∗)–(A7∗), respectively.

  1. (A3∗)

There exist two positive constant $c_{1}<c_{2}$ such that

$c_{1}\leq\lambda_{\min}(I({\bm{B}}))\leq\lambda_{\max}(I({\bm{B}}))\leq
c_{2},$

over the set
$\\{{\bm{B}}:||\mbox{\boldmath$\beta$}_{{\bm{B}}}-\mbox{\boldmath$\beta$}_{{\bm{B}}_{0}}||\leq\triangle\sqrt{p_{n}/n}\\}$
for some constant $\triangle>0$. It is also assumed that $I({\bm{B}})$ has a
constant rank on the same set.

  2. (A4∗)

The true intra-subject correlation matrix ${\bm{R}}_{0}$ has bounded
eigenvalues from zero and infinity. The estimated working correlation matrix
satisfies
$\|\widehat{{\bm{R}}}^{-1}-\tilde{{\bm{R}}}^{-1}\|_{\textbf{F}}=O_{p}(\sqrt{p_{n}/n})$,
where $\|\cdot\|_{\textbf{F}}$ is the Frobenius norm, $\tilde{{\bm{R}}}$ is
some positive definite matrix with bounded eigenvalues from zero and infinity,
and $\tilde{{\bm{R}}}={\bm{R}}_{0}$ is not required.

  3. (A8∗)

It is assumed that $\mu^{(1)}(\theta_{ij})$ are uniformly bounded away from
zero and infinity, and $\mu^{(k)}(\theta_{ij})$ are uniformly bounded by a
finite constant, over the set
$\\{{\bm{B}}:||\mbox{\boldmath$\beta$}_{{\bm{B}}}-\mbox{\boldmath$\beta$}_{{\bm{B}}_{0}}||\leq\triangle\sqrt{p_{n}/n}\\}$,
for some constant $\triangle>0$, $i=1,\dots,n$, $j=1,\dots,m$, and $k=2,3$.

  4. (A9∗)

There exist two positive constants $c_{3}<c_{4}$ such that

$c_{3}\leq\lambda_{\min}({\bm{H}}({\bm{B}},{\bm{X}}_{ij}))\leq\lambda_{\max}({\bm{H}}({\bm{B}},{\bm{X}}_{ij}))\leq
c_{4},$

over the set
$\\{{\bm{B}}:||\mbox{\boldmath$\beta$}_{{\bm{B}}}-\mbox{\boldmath$\beta$}_{{\bm{B}}_{0}}||\leq\triangle\sqrt{p_{n}/n}\\}$
for some constant $\triangle>0$, $i=1,\dots,n$ and $j=1,\dots,m$.

Comparing the two sets of regularity conditions for the fixed and diverging
number of parameters, the main difference is that the conditions are imposed
on the set
$\\{{\bm{B}}:||\mbox{\boldmath$\beta$}_{{\bm{B}}}-\mbox{\boldmath$\beta$}_{{\bm{B}}_{0}}||\leq\triangle\sqrt{p_{n}/n}\\}$
when the number of parameters diverges. This is due to the slower convergence
rate of the tensor GEE estimator with a diverging $p_{n}$. In addition, we
note that $I({\bm{B}})$ and ${\bm{H}}({\bm{B}},{\bm{X}}_{ij})$ are no longer
matrices with fixed dimensions when $p_{n}$ diverges. Correspondingly, we
impose conditions (A3*) and (A9*) on the bounded eigenvalues, which are
similar to the sparse Riesz condition for vector covariates. The latter
condition has been frequently employed in the current literature of inference
with diverging dimensions (Zhang and Huang,, 2008; Zhang,, 2010).

Next we present the asymptotics for the tensor GEE estimator with a diverging
$p_{n}$.

######  Theorem 3.

Under conditions (A1*)-(A9*), and $p_{n}=o(n^{1/2})$, there exists a root
$\widehat{\bm{B}}$ of the equations ${\bm{s}}({\bm{B}})={\mathbf{0}}$
satisfying that

$\displaystyle\|\mbox{\boldmath$\beta$}_{\widehat{{\bm{B}}}}-\mbox{\boldmath$\beta$}_{{\bm{B}}_{0}}\|=O_{p}(\sqrt{p_{n}/n}).$

It is important to note that, if one directly vectorizes the tensor covariate
and applies the asymptotics of the classical GEE as in Wang, (2011), the
conditions for the consistency would require
$\prod_{d=1}^{D}p_{d}=o(n^{1/2})$, i.e. $p_{n}=o(n^{1/(2D)})$. This rate can
be much more stringent for a tensor covariate. Theorem 3, instead, states that
the consistency still holds with $p_{n}=o(n^{1/2})$, after imposing and
exploiting the low rank tensor structure on the coefficients array.

The asymptotic normality can also be established for a diverging $p_{n}$.

######  Theorem 4.

Under conditions (A1*)-(A9*), and $p_{n}=o(n^{1/3})$, for any vector
${\bm{b}}_{n}\in\mathrm{I\\!R}\mathit{{}^{R\sum_{d=1}^{D}p_{d}}}$ such that
$\|{\bm{b}}_{n}\|=1$, we have

$\displaystyle{\bm{b}}_{n}^{\mbox{\tiny{\sf
T}}}\tilde{{\bm{M}}}_{n}^{-1/2}({\bm{B}}_{0})\tilde{{\bm{D}}}_{n1}({\bm{B}}_{0})\bigl{(}\mbox{\boldmath$\beta$}_{\widehat{{\bm{B}}}}-\mbox{\boldmath$\beta$}_{{\bm{B}}_{0}}\bigr{)}\to\mathrm{Normal}(0,1)\textrm{
in distribution.}$

Similarly, for the asymptotic normality to hold, the condition would have
become $p_{n}=o(n^{1/(3D)})$ if one directly vectorizes the tensor covariate.
By contrast, the tensor GEE requires $p_{n}=o(n^{1/3})$.

##  4 Simulations

We have carried out extensive simulations to investigate the finite sample
performance of our proposed tensor GEE approach. We adopt the following
simulation setup. We generated the responses according to the normal model

${\bm{Y}}_{i}\sim\mathrm{MVN}(\mbox{\boldmath$\mu$}_{i},\sigma^{2}{\bf
R}_{0}),\quad i=1,\dots,n,$

where ${\bm{Y}}_{i}=(Y_{i1},\ldots,Y_{im})^{\mbox{\tiny{\sf T}}}$,
$\mbox{\boldmath$\mu$}_{i}=(\mu_{i1},\ldots,\mu_{im})^{\mbox{\tiny{\sf T}}}$,
$\sigma^{2}$ is a scale parameter, and ${\bm{R}}_{0}$ is the true $m\times m$
intra-subject correlation matrix. We have chosen ${\bm{R}}_{0}$ to be of an
exchangeable (compound symmetric) structure with the off-diagonal coefficient
$\rho=0.8$. The mean function is of the form

$\mu_{ij}=\mbox{\boldmath$\gamma$}^{\mbox{\tiny{\sf
T}}}{\bm{Z}}_{ij}+\langle{\bm{B}},{\bm{X}}_{ij}\rangle,\quad
i=1,\ldots,n,j=1,\ldots,m,$

where ${\bm{Z}}_{ij}\in\mathrm{I\\!R}\mathit{{}^{5}}$ denotes the covariate
vector, with all elements generated from a standard normal distribution, and
$\mbox{\boldmath$\gamma$}\in\mathrm{I\\!R}\mathit{{}^{5}}$ is the
corresponding coefficient vector, with all elements equal to one;
${\bm{X}}_{ij}\in\mathrm{I\\!R}\mathit{{}^{64\times 64}}$ denotes the 2D
matrix covariate, again with all elements from standard normal, and
${\bm{B}}\in\mathrm{I\\!R}\mathit{{}^{64\times 64}}$ is the matrix
coefficient. ${\bm{B}}$ takes the value of 0 or 1, and contains a series of
shapes as shown in Figure 1, including “square”, “T-shape”, “disk”,
“triangle”, and “butterfly”. Our goal is to recover those shapes in ${\bm{B}}$
by inferring the association between $Y_{ij}$ and ${\bm{X}}_{ij}$ after
adjusting for ${\bm{Z}}_{ij}$.

###  4.1 Signal Recovery

As the true signal in reality is hardly of an exact low rank structure, the
tensor model (4) and the associated tensor GEE (5) essentially provide a low
rank _approximation_ to the true signal. It is thus important to verify if
such an approximation is adequate. We set $n=500$, $m=4$, and show both the
tensor GEE estimates under various ranks and the corresponding BIC values (7)
in Figure 1. We first assume that the correlation structure is correctly
specified, and will study potential misspecification in the next section. In
this setup, “square” has the true rank equal to 1, “T-shape” has the rank 2,
and the remaining shapes have the highest possible rank 64. It is clearly seen
from the figure that the tensor GEE offers a sound recovery of the true
signal, even for the signals with high rank or natural shape, e.g., “disk” and
“butterfly”. In addition, the BIC seems to identify the correct or best
approximate rank for all the signals.

---
Figure 1: True and recovered image signals by the tensor GEE with varying
ranks. $n=500,m=4$. The correlation structure is correctly specified. TR$(R)$
means estimate from the rank-$R$ tensor model.

###  4.2 Effect of Correlation Specification

Table 1: Bias, variance, and MSE of the tensor GEE estimates under various working correlation structures. Reported are the average out of 100 simulation replicates. The true intra-subject correlation is exchangeable with $\rho=0.8$. $n$ | $m$ | Working Correlation | Bias2 | Variance | MSE
---|---|---|---|---|---
50 | 10 | Exchangeable | 122.0 | 383.6 | 505.6(7.9)
|  | AR-1 | 139.1 | 530.0 | 669.1(15.8)
|  | Independence | 119.1 | 393.9 | 513.0(11.0)
100 | 10 | Exchangeable | 85.8 | 128.9 | 214.7(2.2)
|  | AR-1 | 88.0 | 159.1 | 247.1(3.0)
|  | Independence | 93.0 | 141.2 | 234.2(2.8)
150 | 10 | Exchangeable | 86.1 | 51.3 | 137.2(0.6)
|  | AR-1 | 85.6 | 56.0 | 141.6(0.6)
|  | Independence | 84.9 | 62.3 | 147.2(0.9)

We have shown that the tensor GEE estimator remains asymptotically consistent
even when the working correlation structure is misspecified. However this
describes only the _large sample_ behavior. In this section, we investigate
potential effect of correlation misspecification when the sample size is
_small_ or _moderate_.

We chose the “butterfly” signal and fitted the tensor GEE model with three
different working correlation structures: exchangeable, which is the correct
specification in our setup, autoregressive of order one (AR-1), and
independent. Table 1 reports the averages and standard errors out of 100
replicates of the squared bias, the variance, and the mean squared error (MSE)
of the tensor GEE estimate. We observe that the estimator based on the correct
working correlation structure, i.e., the exchangeable structure, performs
better than those based on misspecified correlation structures. When the
sample size is moderate ($n=100$), all the estimators have comparable bias,
while the difference in MSE mostly comes from the variance part of the
estimator. This agrees with the theory that the choice of the working
correlation structure affects the asymptotic variance of the estimator. When
the sample size becomes relatively large ($n=150$), all the estimators perform
similarly by the scaling term of $n^{-1/2}$ on the variance. When the sample
size is small ($n=50$), all the estimators have relatively large bias, while
the independence working structure yield similar results as the exchangeable
structure. This suggests that, when the sample size is _limited_ , using a
simple independence working structure is probably preferable compared to a
more complex correlation structure.

---
Figure 2: Snapshots of tensor GEE estimation with different working
correlation structures. The true correlation is an equicorrelated structure.
The comparison is _row-wise_. The first row shows a replicate where the
estimates are “close” to the average behavior, and thus the visual quality of
the estimates under different correlations structures are similar. The second
row shows a replicate where the estimates are “far away” from the average,
then the estimate under the correct correlation structure (panel 1) is clearly
superior than those under incorrect structures.

Nevertheless, we should bear in mind that the above observations are for the
_average_ behavior of the estimate. Figure 2 shows _two snapshots_ of the
estimated signals under the three working correlations at $n=100$. The top
panel is one replicate where the estimates are “close” to the average in the
sense that the bias, variance and MSE values for this single data realization
are similar to those averages reported in Table 1. Consequently, the visual
qualities of the three recovered signals are similar. The bottom panel, on the
other hand, shows another replicate where the estimates are “far away” from
the average. Then the quality of the estimated signal under the correct
working correlation structure is superior than the ones under the incorrect
specifications. Such an observation suggests that, as long as the sample size
of the study is moderate to large, a longitudinal model should be favored over
the one that totally ignores potential intra-subject correlation.

###  4.3 Regularized Estimation

We implemented the regularized tensor GEE with a lasso penalty, which extends
the penalized GEE method of Wang et al., (2012) from vector to array
covariate. It can identify relevant regions in images that are associated with
the outcome, and this _region selection_ problem corresponds to _variable
selection_ in classical vector covariate regressions. We studied the empirical
performance by adopting the simulation setup described at the beginning of
Section 4, but varying the sample size. The estimates of three shapes,
“T-shape”, “triangle”, and “butterfly”, with and without regularizations, are
shown in Figure 3. For the regularized tensor GEE, the penalty parameter
$\lambda$ was selected based on the prediction accuracy on an independent
validation set. It is clearly seen from the plot that, while increasing sample
size improves estimation accuracy for both tensor GEE and regularized tensor
GEE, regularization leads to a more accurate recovery, especially when the
sample size is limited. As such we recommend the regularized tensor GEE for
longitudinal imaging data analysis in practice.

---
Figure 3: Comparison of tensor GEE estimation with and without regularization
under varying sample size. $m=4$. The matrix covariate is of size $64\times
64$.

##  5 Real Data Analysis

###  5.1 Alzheimer’s Disease

Alzheimer’s Disease (AD) is a progressive and irreversible neurodegenerative
disorder and the leading form of dementia in elderly subjects. It is
characterized by gradual impairment of cognitive and memory functions, and it
has been projected to quadruple in its prevalence by the year 2050 (Brookmeyer
et al.,, 2007). Amnestic mild cognitive impairment (MCI) is often a prodromal
stage to Alzheimer’s disease, and individuals with MCI may convert to AD at an
annual rate as high as $15\%$ (Petersen et al.,, 1999). As such there is a
pressing need for accurate and early diagnosis of AD and MCI, as well as
monitoring their progression. The data we analyzed was obtained from the
Alzheimer’s Disease Neuroimaging Initiative (ADNI). It consists of $n=88$ MCI
subjects with longitudinal MRI images of white matter at baseline, 6-month,
12-month, 18-month and 24-month ($m=5$). Also recorded for each subject at
multiple visits was the Mini Mental State Examination (MMSE) score. It
measures the orientation to time and place, the immediate and delayed recall
of three words, the attention and calculations, language, and
visuoconstructional functions (Folstein et al.,, 1975), and is our response
variable. A detailed description of acquiring MRI data from ADNI and the
preprocessing protocol can be found in Zhang et al., (2012). There are two
scientific goals for this study. One is to predict the future clinical scores
based on the data at previous time points, which is particularly useful for
monitoring disease progression. The second is to identify brain subregions
that are highly relevant to the disorder. We fitted tensor GEE to this data
for both score prediction and region selection.

###  5.2 Prediction and Disease Prognosis

We downsized the original $256\times 256\times 256$ MRI images to $32\times
32\times 32$ via interpolation for computational simplicity. We first fitted
tensor GEE using the data from baseline to 12-month, and used prediction of
MMSE at 18-month to select the tuning parameter $\lambda$. Then we refitted
the model using the data from baseline to 18-month under the selected
$\lambda$, and evaluated the prediction accuracy of all subjects using the
“future” MMSE score at 24-month. The accuracy was evaluated by the rooted mean
squared error (RMSE),
$\\{n^{-1}\sum_{i=1}^{n}(Y_{im}-\hat{Y}_{im})^{2}\\}^{1/2}$, and the
correlation, $\textrm{Corr}(Y_{im},\hat{Y}_{im})$. This evaluation scheme is
the same as that of Zhang et al., (2012). Table 2 summarizes the results. It
is seen that, for this data set, the best prediction was achieved under an
AR(1) working correlation structure with $L_{1}$ regularization. The
corresponding RMSE and correlation were 2.270 and 0.747, which are only
slightly worse than the best reported RMSE 2.035 and correlation 0.786 in
Zhang et al., (2012). Note that Zhang et al., (2012) used multiple imaging
modalities and additional clinical covariates, which are supposed to improve
the prediction accuracy, while our study utilized only one imaging modality.

Table 2: Prediction of future clinical MMSE scores using tensor GEE | RMSE: $\\{\sum_{i=1}^{n}n^{-1}(Y_{im}-\hat{Y}_{im})^{2}\\}^{1/2}$
---|---
Working Correlation | Independence | Equicorrelated | AR(1) | Unstructured
regularization | 2.460 | 2.349 | 2.270 | 2.570
no regularization | 2.526 | 2.427 | 2.429 | 2.628
| Correlation: $\textrm{Corr}(Y_{im},\hat{Y}_{im})$
Working Correlation | Independence | Equicorrelated | AR(1) | Unstructured
regularization | 0.705 | 0.733 | 0.747 | 0.700
no regularization | 0.701 | 0.716 | 0.725 | 0.693

###  5.3 Region Selection

We applied the lasso regularized tensor GEE to this data, and Figure 4 shows
the estimate (marked in red) overlaid on an image of an arbitrarily chosen
subject, with three views, top, side and bottom, respectively. The identified
anatomical regions mainly correspond to cerebral cortex, part of temporal
lobe, parietal lobe, and frontal lobe (Braak and Braak,, 1991; Desikan et
al.,, 2009; Yao et al.,, 2012). With AD, patients experience significant
widespread damage over the brain, causing shrinkage of brain volume (Yao et
al.,, 2012; Harasty et al.,, 1999) and thinning of cortical thickness (Desikan
et al.,, 2009; Yao et al.,, 2012). The affected brain regions include those
involved in controlling language (Broca’s area) (Harasty et al.,, 1999),
reasoning (superior and inferior frontal gyri) (Harasty et al.,, 1999), part
of sensory area (primary auditory cortex, olfactory cortex, insula, and
operculum) (Braak and Braak,, 1991; Lee et al.,, 2013), somatosensory
association area (Yao et al.,, 2012; Tales et al.,, 2005; Mapstone et al.,,
2003), memory loss (hippocampus) (den Heijer et al.,, 2010), and motor
function (Buchman and Bennett,, 2011). It is interesting to note that these
regions are affected starting at different stages of AD, indicating the
capability of the proposed method to locate brain atrophies as disease
progresses. Specifically, hippocampus, which is highly correlated to memory
loss, is commonly detected at the earliest stage of the disease. Regions
related to language, communication, and motor functions are normally detected
at the later stages of the disease. The fact that our findings are consistent
with the results reported in previous studies demonstrates the efficacy of our
proposed method in identifying correct biomarkers that are closely related to
AD/MCI.

|  |
---|---|---
Figure 4: The ADNI data: regularized estimate overlaid on a randomly selected
subject.

##  6 Discussions

We have proposed a tensor GEE approach for analyzing longitudinal imaging
data. Our method combines the powerful GEE idea for handling longitudinal
correlation and the low rank tensor decomposition to reduce the vast
dimensionality of imaging data. The proposed algorithm scales well with
imaging data size and is easy to implement using existing statistical
softwares. Simulation studies and real data analysis show the advantage of our
method for both signal recovering and prediction.

In the current paper, we have considered an image covariate together with a
conventional vector covariate. Extending to joint multi-modality imaging
analysis is conceptually easy: simply adding more array covariates into the
systematic component (4). However this brings up other issues such as joint
selection of ranks for multiple array covariates, properly defining
interactions between tensor covariates, and even higher volume of data. These
important yet nontrivial questions deserve further investigation.

## References

  * ADHD, (2014) ADHD (2014).  The ADHD-200 sample.  http://fcon_1000.projects.nitrc.org/indi/adhd200/.  [Online; accessed 04-Mar-2014].
  * ADNI, (2014) ADNI (2014).  Alzheimer’s disease neuroimaging initiative.  http://adni.loni.usc.edu.  [Online; accessed 04-Mar-2014].
  * Aston and Kirch, (2012) Aston, J. A. and Kirch, C. (2012).  Estimation of the distribution of change-points with application to fmri data.  Annals of Applied Statistics, 6:1906–1948.
  * Balan and Schiopu-Kratina, (2005) Balan, R. M. and Schiopu-Kratina, I. (2005).  Asymptotic results with generalized estimating equations for longitudinal data.  The Annals of Statistics, 33(2):522–541.
  * Braak and Braak, (1991) Braak, H. and Braak, E. (1991).  Neuropathological stageing of Alzheimer-related changes.  Acta Neuropathologica, 82(4):239–259.
  * Brookmeyer et al., (2007) Brookmeyer, R., Johnson, E., Ziegler-Graham, K., and Arrighi, H. M. (2007).  Forecasting the global burden of alzheimer’s disease.  Alzheimer’s & Dementia, 3(3):186 – 191.
  * Buchman and Bennett, (2011) Buchman, A. and Bennett, D. (2011).  Loss of motor function in preclinical alzheimer’s disease.  Expert Review Neurotherapeutics, 11(5):665–676.
  * Caffo et al., (2010) Caffo, B., Crainiceanu, C., Verduzco, G., Joel, S., S.H., M., Bassett, S., and Pekar, J. (2010).  Two-stage decompositions for the analysis of functional connectivity for fMRI with application to Alzheimer’s disease risk.  NeuroImage, 51(3):1140–1149.
  * Davatzikos et al., (2009) Davatzikos, C., Xu, F., An, Y., Fan, Y., and Resnick, S. M. (2009).  Longitudinal progression of alzheimer’s-like patterns of atrophy in normal older adults: the spare-ad index.  Brain, 132(8):2026–2035.
  * den Heijer et al., (2010) den Heijer, T., van der Lijn, F., Koudstaal, P. J., Hofman, A., van der Lugt, A., Krestin, G. P., Niessen, W. J., and Breteler, M. M. B. (2010).  A 10-year follow-up of hippocampal volume on magnetic resonance imaging in early dementia and cognitive decline.  Brain, 133(4):1163–1172.
  * Desikan et al., (2009) Desikan, R., Cabral, H., Hess, C., Dillon, W., Salat, D., Buckner, R., Fischl, B., and Initiative, A. D. N. (2009).  Automated MRI measures identify individuals with mild cognitive impairment and Alzheimer’s disease.  Brain, 132:2048–2057.
  * Fan and Li, (2001) Fan, J. and Li, R. (2001).  Variable selection via nonconcave penalized likelihood and its oracle properties.  J. Amer. Statist. Assoc., 96(456):1348–1360.
  * Fan and Li, (2004) Fan, J. and Li, R. (2004).  New estimation and model selection procedures for semiparametric modeling in longitudinal data analysis.  Journal of the American Statistical Association, 99(467):710–723.
  * Folstein et al., (1975) Folstein, M. F., Folstein, S. E., and McHugh, P. R. (1975).  Mini-mental state: A practical method for grading the cognitive state of patients for the clinician.  Journal of Psychiatric Research, 12(3):189 – 198.
  * Frank and Friedman, (1993) Frank, I. E. and Friedman, J. H. (1993).  A statistical view of some chemometrics regression tools.  Technometrics, 35(2):109–135.
  * Friston, (2009) Friston, K. J. (2009).  Modalities, modes, and models in functional neuroimaging.  Science, 326:399–403.
  * Harasty et al., (1999) Harasty, J. A., Halliday, G. M., Kril, J. J., and Code, C. (1999).  Specific temporoparietal gyral atrophy reflects the pattern of language dissolution in alzheimer’s disease.  Brain, 122(4):675–686.
  * Hinrichs et al., (2009) Hinrichs, C., Singh, V., Mukherjee, L., Xu, G., Chung, M. K., Johnson, S. C., and ADNI (2009).  Spatially augmented lpboosting for ad classification with evaluations on the adni dataset.  NeuroImage, 48:138–149.
  * Hinrichs et al., (2011) Hinrichs, C., Singh, V., Xu, G., and Johnson, S. C. (2011).  Predictive markers for {AD} in a multi-modality framework: An analysis of {MCI} progression in the {ADNI} population.  NeuroImage, 55(2):574 – 589.
  * Kang et al., (2012) Kang, H., Ombao, H., Linkletter, C., Long, N., and Badre, D. (2012).  Spatio-spectral mixed-effects model for functional magnetic resonance imaging data.  Journal of the American Statistical Association, 107(498):568–577.
  * Kolda and Bader, (2009) Kolda, T. G. and Bader, B. W. (2009).  Tensor decompositions and applications.  SIAM Rev., 51(3):455–500.
  * Lazar, (2008) Lazar, N. A. (2008).  The Statistical Analysis of Functional MRI Data.  Springer, New York.
  * Lee et al., (2013) Lee, T. M., Sun, D., Leung, M.-K., Chu, L.-W., and Keysers, C. (2013).  Neural activities during affective processing in people with alzheimer’s disease.  Neurobiology of Aging, 34(3):706 – 715.
  * Li, (1997) Li, B. (1997).  On the consistency of generalized estimating equations.  In Selected Proceedings of the Symposium on Estimating Functions (Athens, GA, 1996), volume 32 of IMS Lecture Notes Monogr. Ser., pages 115–136. Inst. Math. Statist., Hayward, CA.
  * Li et al., (2013) Li, Y., Gilmore, J. H., Shen, D., Styner, M., Lin, W., and Zhu, H. (2013).  Multiscale adaptive generalized estimating equations for longitudinal neuroimaging data.  NeuroImage, 72(0):91 – 105.
  * Liang and Zeger, (1986) Liang, K. Y. and Zeger, S. L. (1986).  Longitudinal data analysis using generalized linear models.  Biometrika, 73(1):13–22.
  * Mapstone et al., (2003) Mapstone, M., Steffenella, T., and Duffy, C. (2003).  A visuospatial variant of mild cognitive impairment: getting lost between aging and ad.  Neurology, 60:802–808.
  * McEvoy et al., (2011) McEvoy, L. K., Holland, D., Hagler, D. J., Fennema-Notestine, C., Brewer, J. B., and Dale, A. M. (2011).  Mild cognitive impairment: Baseline and longitudinal structural mr imaging measures improve predictive prognosis.  Radiology, 259(3):834–843.  PMID: 21471273.
  * Misra et al., (2009) Misra, C., Fan, Y., and Davatzikos, C. (2009).  Baseline and longitudinal patterns of brain atrophy in {MCI} patients, and their use in prediction of short-term conversion to ad: Results from {ADNI}.  NeuroImage, 44(4):1415 – 1422.
  * Ni et al., (2010) Ni, X., Zhang, D., and Zhang, H. H. (2010).  Variable selection for semiparametric mixed models in longitudinal studies.  Biometrics, 66(1):79–88.
  * Ortega and Rheinboldt, (2000) Ortega, J. M. and Rheinboldt, W. C. (2000).  Iterative solution of nonlinear equations in several variables, volume 30.  Siam.
  * Pan, (2001) Pan, W. (2001).  Akaike’s information criterion in generalized estimating equations.  Biometrics, 57(1):120–125.
  * Petersen et al., (1999) Petersen, R., Smith, G., Waring, S., Ivnik, R., Tangalos, E., and Kokmen, E. (1999).  Mild cognitive impairment: clinical characterization and outcome.  Archives of Neurology, 56:303–308.
  * Prentice and Zhao, (1991) Prentice, R. L. and Zhao, L. P. (1991).  Estimating equations for parameters in means and covariances of multivariate discrete and continuous responses.  Biometrics, 47(3):825–839.
  * Qu et al., (2000) Qu, A., Lindsay, B. G., and Li, B. (2000).  Improving generalised estimating equations using quadratic inference functions.  Biometrika, 87(4):823–836.
  * Rao and Mitra, (1971) Rao, C. R. and Mitra, S. K. (1971).  Generalized Inverse of Matrices and its Applications.  John Wiley & Sons, Inc., New York-London-Sydney.
  * Reiss and Ogden, (2010) Reiss, P. and Ogden, R. (2010).  Functional generalized linear models with images as predictors.  Biometrics, 66:61–69.
  * Skup et al., (2012) Skup, M., Zhu, H., and Zhang, H. (2012).  Multiscale adaptive marginal analysis of longitudinal neuroimaging data with time-varying covariates.  Biometrics, 68(4):1083–1092.
  * Song et al., (2009) Song, P. X.-K., Jiang, Z., Park, E., and Qu, A. (2009).  Quadratic inference functions in marginal models for longitudinal data.  Statistics in Medicine, 28(29):3683–3696.
  * Tales et al., (2005) Tales, A., Haworth, J., Nelson, S., J. Snowden, R., and Wilcock, G. (2005).  Abnormal visual search in mild cognitive impairment and alzheimer’s disease.  Neurocase, 11(1):80–84.
  * Tibshirani, (1996) Tibshirani, R. (1996).  Regression shrinkage and selection via the lasso.  J. Roy. Statist. Soc. Ser. B, 58(1):267–288.
  * Wang, (2011) Wang, L. (2011).  GEE analysis of clustered binary data with diverging number of covariates.  The Annals of Statistics, 39(1):389–417.
  * Wang et al., (2012) Wang, L., Zhou, J., and Qu, A. (2012).  Penalized generalized estimating equations for high-dimensional longitudinal data analysis.  Biometrics, 68(2):353–360.
  * Wang et al., (2014) Wang, X., Nan, B., Zhu, J., and Koeppe, R. (2014).  Regularized 3D functional regression for brain image data via haar wavelets.  The Annals of Applied Statistics, page in press.
  * Xie and Yang, (2003) Xie, M. and Yang, Y. (2003).  Asymptotics for generalized estimating equations with large cluster sizes.  The Annals of Statistics, 31(1):310–347.
  * Xue et al., (2010) Xue, L., Qu, A., and Zhou, J. (2010).  Consistent model selection for marginal generalized additive model for correlated data.  Journal of the American Statistical Association, 105(492):1518–1530.  Supplementary materials available online.
  * Yao et al., (2012) Yao, Z., Hu, B., Liang, C., Zhao, L., Jackson, M., and the Alzheimer’s Disease Neuroimaging Initiative (2012).  A longitudinal study of atrophy in amnestic mild cognitive impairment and normal aging revealed by cortical thickness.  PLoS One, 7(11):e48973.
  * Zhang, (2010) Zhang, C.-H. (2010).  Nearly unbiased variable selection under minimax concave penalty.  The Annals of Statistics, 38(2):894–942.
  * Zhang and Huang, (2008) Zhang, C.-H. and Huang, J. (2008).  The sparsity and bias of the lasso selection in high-dimensional linear regression.  The Annals of Statistics, pages 1567–1594.
  * Zhang et al., (2012) Zhang, D., Shen, D., and Alzheimer’s Disease Neuroimaging Initiative (2012).  Predicting future clinical changes of mci patients using longitudinal and multimodal biomarkers.  PLoS One, 7(3):e33182.
  * Zhang et al., (2011) Zhang, D., Wang, Y., Zhou, L., Yuan, H., Shen, D., and the Alzheimers Disease Neuroimaging Initiative (2011).  Multimodal classification of Alzheimer’s disease and mild cognitive impairment.  NeuroImage, 55(3):856 – 867.
  * Zhou et al., (2013) Zhou, H., Li, L., and Zhu, H. (2013).  Tensor regression with applications in neuroimaging data analysis.  Journal of the American Statistical Association, 108(502):540–552.
  * Zou and Hastie, (2005) Zou, H. and Hastie, T. (2005).  Regularization and variable selection via the elastic net.  J. R. Stat. Soc. Ser. B Stat. Methodol., 67(2):301–320.

## Appendix: Technical Proofs

Outline of the proofs
We prove the results for the diverging case (Theorem 3 and Theorem 4) in the
appendix. One can prove the results for the fixed case (Theorem 1 and Theorem
2) by using the same techniques below and replacing $p_{n}$ with a fixed
positive constant.

The proof of Lemma 1 is similar to the one of Theorem 3 by dropping the terms
involving the working correlation matrix and thus is omitted here.

To facilitate the proof, we introduce the following notations. Denote
$\widehat{\bm{\beta}}_{n}=\bm{\beta}_{\hat{{\bm{B}}}}$ and
$\bm{\beta}_{0}=\mbox{\boldmath$\beta$}_{{\bm{B}}_{0}}$. Recall that the CP
decomposition ensures that ${\bm{B}}$ is uniquely determined by
$\bm{\beta}_{n}\in\mathrm{I\\!R}\mathit{{}^{R\sum_{d=1}^{D}p_{d}}}$. Denote
${\bm{J}}(\bm{\beta})=[{\bm{J}}_{1}\cdots{\bm{J}}_{D}]$, and note that under
tensor structure
$\partial\theta_{ij}/\partial\bm{\beta}={\bm{J}}(\bm{\beta})\text{vec}({\bm{X}}_{ij})$.
Recall the generalized estimating equations without vector covariates can be
written as

$\displaystyle{\bm{s}}_{n}(\bm{\beta}_{n})$
$\displaystyle=\sum_{i=1}^{n}{\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{n})\text{vec}({\bm{X}}_{i}){\bm{A}}^{1/2}_{i}(\bm{\beta}_{n})\widehat{{\bm{R}}}^{-1}{\bm{A}}^{-1/2}_{i}(\bm{\beta}_{n})({\bm{Y}}_{i}-\bm{\mu}_{i}(\bm{\beta}_{n})).$

The main technique to prove Theorem 3 is the sufficient condition for
existence and consistency of a root of equations proposed in Ortega and
Rheinboldt, (2000). To check this condition, the following Lemma 2 \- 4 are
proposed. Lemma 2 provides a useful approximation to the generalized
estimating equations ${\bm{s}}_{n}(\bm{\beta}_{0})$ based on the Condition
(A4*) of the working correlation matrix. This facilitates the later
evaluations of the moments of the generalized estimating equations by treating
the intra-subject correlation as known. Lemma 3 further establishes the
approximation to the negative gradients of the generalized estimating
equations. Lemma 4 refines this approximation to the negative gradients at one
more step, providing the foundations for the Talyor expansion of generalized
estimating equations at the true value.

Based on Theorem 3, the proof of Theorem 4 is straightforward by evaluating
the covariance matrix of the generalized estimating equations and applying the
Lindeberg-Feller central limit theorem.

######  Lemma 2.

Under Conditions (A1*)-(A9*), $p_{n}=o(n^{1/2})$, then
$||\tilde{{\bm{s}}}_{n}(\bm{\beta}_{0})-{\bm{s}}_{n}(\bm{\beta}_{0})||=O_{p}(p_{n})$,
where $\tilde{{\bm{s}}}_{n}(\bm{\beta}_{0})$ is ${\bm{s}}_{n}(\bm{\beta}_{0})$
with $\widehat{{\bm{R}}}$ replaced by $\tilde{{\bm{R}}}$.

###### Proof of Lemma 2.

Consider

$\displaystyle\tilde{{\bm{s}}}_{n}(\bm{\beta}_{n})$
$\displaystyle=\sum_{i=1}^{n}{\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{n})\text{vec}({\bm{X}}_{i}){\bm{A}}^{1/2}_{i}(\bm{\beta}_{n})\tilde{{\bm{R}}}^{-1}{\bm{A}}^{-1/2}_{i}(\bm{\beta}_{n})({\bm{Y}}_{i}-\bm{\mu}_{i}(\bm{\beta}_{n})).$

Denote $\\{r_{i,j}\\}_{1\leq i,j\leq m}$ the $(i,j)$-th element of
$\widehat{{\bm{R}}}^{-1}-\tilde{{\bm{R}}}^{-1}$. By Condition (A4*),
$r_{i,j}=O_{p}(\sqrt{p_{n}/n})$. Note that

$\displaystyle{\bm{s}}_{n}(\bm{\beta}_{0})-\tilde{{\bm{s}}}_{n}(\bm{\beta}_{0})$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{n}\sum_{j=1}^{m}\sum_{k=1}^{m}r_{j,m}\sigma_{ij}(\bm{\beta}_{0})\epsilon_{ik}(\bm{\beta}_{0}){\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{0})\text{vec}({\bm{X}}_{ij})$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{m}\sum_{k=1}^{m}r_{j,m}\Big{[}\sum_{i=1}^{n}\sigma_{ij}(\bm{\beta}_{0})\epsilon_{ik}(\bm{\beta}_{0}){\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{0})\text{vec}({\bm{X}}_{ij})\Big{]},$

where
$\epsilon_{ik}(\bm{\beta}_{0})=\sigma_{ik}^{-1}(\bm{\beta}_{0})(Y_{ik}-\mu_{ik}(\bm{\beta}_{0}))$.
By Condition (A6*), $\mathbb{E}[\epsilon_{ik}(\bm{\beta}_{0})]=O_{p}(1)$. Note
that for any $1\leq j,k\leq m$,

$\displaystyle\mathbb{E}\Big{[}||\sum_{i=1}^{n}\sigma_{ij}(\bm{\beta}_{0})\epsilon_{ik}(\bm{\beta}_{0}){\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{n})\text{vec}({\bm{X}}_{ij})||^{2}\Big{]}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{n}\sigma^{2}_{ij}(\bm{\beta}_{0})\mathbb{E}[\epsilon_{ik}^{2}(\bm{\beta}_{0})]\text{vec}^{\mbox{\tiny{\sf
T}}}({\bm{X}}_{ij}){\bm{J}}(\bm{\beta}_{0}){\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{0})\text{vec}({\bm{X}}_{ij})$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{n}\sigma^{2}_{ij}(\bm{\beta}_{0})\mathbb{E}[\epsilon_{ik}^{2}(\bm{\beta}_{0})]\text{Tr}({\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{0})\text{vec}({\bm{X}}_{ij})\text{vec}^{\mbox{\tiny{\sf
T}}}({\bm{X}}_{ij}){\bm{J}}(\bm{\beta}_{0}))$ $\displaystyle\leq$
$\displaystyle Cnp_{n},$

for some constant $C>0$ by Condition (A1*), (A2*) and (A7*). Since
$r_{i,j}=O_{p}(\sqrt{p_{n}/n})$, the proof is complete. ∎

Consider
${\bm{D}}_{n}(\bm{\beta}_{n})=-\partial{\bm{s}}_{n}(\bm{\beta}_{n})/\partial\bm{\beta}_{n}$,
$\tilde{{\bm{D}}}_{n}(\bm{\beta}_{n})=-\partial\tilde{{\bm{s}}}_{n}(\bm{\beta}_{n})/\partial\bm{\beta}_{n}$.
Lemma 3 establishes the approximation to the negative gradients of the
estimating equations.

######  Lemma 3.

Under Conditions (A1*)-(A9*), for any $\triangle>0$,

$\displaystyle\sup_{||\bm{\beta}_{n}-\bm{\beta}_{0}||\leq\triangle\sqrt{p_{n}/n}}|\lambda_{\max}[\tilde{{\bm{D}}}_{n}(\bm{\beta}_{n})-{\bm{D}}_{n}(\bm{\beta}_{n})]|=O_{p}(\sqrt{p_{n}n}),$
$\displaystyle\sup_{||\bm{\beta}_{n}-\bm{\beta}_{0}||\leq\triangle\sqrt{p_{n}/n}}|\lambda_{\min}[\tilde{{\bm{D}}}_{n}(\bm{\beta}_{n})-{\bm{D}}_{n}(\bm{\beta}_{n})]|=O_{p}(\sqrt{p_{n}n}).$

###### Proof of Lemma 3.

Similar to Lemma C.1. of Wang, (2011), it can be shown by direct calculation
that

$\tilde{{\bm{D}}}_{n}(\bm{\beta}_{n})=\tilde{{\bm{D}}}_{n1}(\bm{\beta}_{n})+\tilde{{\bm{D}}}_{n2}(\bm{\beta}_{n})+\tilde{{\bm{D}}}_{n3}(\bm{\beta}_{n})+\tilde{{\bm{D}}}_{n4}(\bm{\beta}_{n}),$

where

$\displaystyle\tilde{{\bm{D}}}_{n1}(\bm{\beta}_{n})$
$\displaystyle=\sum_{i=1}^{n}{\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{n})\text{vec}({\bm{X}}_{i}){\bm{A}}^{1/2}_{i}(\bm{\beta}_{n})\tilde{{\bm{R}}}^{-1}{\bm{A}}^{1/2}_{i}(\bm{\beta}_{n})\text{vec}^{\mbox{\tiny{\sf
T}}}({\bm{X}}_{i}){\bm{J}}(\bm{\beta}_{n}),$
$\displaystyle\tilde{{\bm{D}}}_{n2}(\bm{\beta}_{n})$
$\displaystyle=\frac{1}{2}\sum_{i=1}^{n}{\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{n})\text{vec}({\bm{X}}_{i}){\bm{A}}^{1/2}_{i}(\bm{\beta}_{n})\tilde{{\bm{R}}}^{-1}{\bm{A}}^{-3/2}_{i}(\bm{\beta}_{n}){\bf
C}_{i}(\bm{\beta}_{n}){\bm{F}}_{i}(\bm{\beta}_{n})\text{vec}^{\mbox{\tiny{\sf
T}}}({\bm{X}}_{i}){\bm{J}}(\bm{\beta}_{n}),$
$\displaystyle\tilde{{\bm{D}}}_{n3}(\bm{\beta}_{n})$
$\displaystyle=-\frac{1}{2}\sum_{i=1}^{n}{\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{n})\text{vec}({\bm{X}}_{i}){\bm{A}}^{1/2}_{i}(\bm{\beta}_{n}){\bm{F}}_{i}(\bm{\beta}_{n}){\bm{K}}_{i}(\bm{\beta}_{n})\text{vec}^{\mbox{\tiny{\sf
T}}}({\bm{X}}_{i}){\bm{J}}(\bm{\beta}_{n}),$
$\displaystyle\tilde{{\bm{D}}}_{n4}(\bm{\beta}_{n})$
$\displaystyle=\sum_{i=1}^{n}\sum_{j=1}^{m}{\bf e}_{j}^{\mbox{\tiny{\sf
T}}}{\bm{A}}^{1/2}_{i}(\bm{\beta}_{n})\tilde{{\bm{R}}}^{-1}{\bm{A}}^{-1/2}_{i}(\bm{\beta}_{n})({\bm{Y}}_{i}-\bm{\mu}_{i}(\bm{\beta}_{n})){\bf
H}(\bm{\beta}_{n},{\bm{X}}_{ij}),$

with

$\displaystyle{\bf C}_{i}(\bm{\beta}_{n})$
$\displaystyle=\text{diag}\Big{(}Y_{i1}-\mu_{i1}(\bm{\beta}_{n}),\dots,Y_{im}-\mu_{im}(\bm{\beta}_{n})\Big{)},$
$\displaystyle{\bm{F}}_{i}(\bm{\beta}_{n})$
$\displaystyle=\text{diag}\Big{(}\mu_{i1}^{(2)}(\bm{\beta}_{n}),\dots,\mu_{im}^{(2)}(\bm{\beta}_{n})\Big{)},$
$\displaystyle{\bm{K}}_{i}(\bm{\beta}_{n})$
$\displaystyle=\text{diag}\Big{(}\tilde{{\bm{R}}}^{-1}{\bm{A}}^{-1/2}_{i}(\bm{\beta}_{n})({\bm{Y}}_{i}-\bm{\mu}_{i}(\bm{\beta}_{n}))\Big{)},$

${\bf e}^{\mbox{\tiny{\sf T}}}_{j}$ length $m$ vector with $j$-th element 1
and 0 everywhere else, and ${\bf
H}(\bm{\beta}_{n},{\bm{X}}_{ij})=\partial{\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{n})\text{vec}({\bm{X}}_{ij})/\partial\bm{\beta}^{\mbox{\tiny{\sf
T}}}_{n}$.

Let ${\bm{D}}_{ni}(\bm{\beta}_{n})$ be defined the same as
$\tilde{{\bm{D}}}_{ni}(\bm{\beta}_{n})$, but with $\tilde{{\bm{R}}}$ replaced
by $\widehat{{\bm{R}}}$, for $i=1,\dots,4$. It is sufficient to prove

$\sup_{||\bm{\beta}_{n}-\bm{\beta}_{0}||\leq\triangle\sqrt{p_{n}/n}}\sup_{{\bm{u}}}|{\bm{u}}^{\mbox{\tiny{\sf
T}}}[{\bm{D}}_{ni}(\bm{\beta}_{n})-\tilde{{\bm{D}}}_{ni}(\bm{\beta}_{n})]{\bm{u}}|=O_{p}(\sqrt{p_{n}n})$

for any ${\bm{u}}\in\mathrm{I\\!R}\mathit{{}^{R\sum_{d=1}^{D}p_{d}}}$ such
that $||{\bm{u}}||=1$, $i=1,\dots,4$.

For $i=1$, we have

$\displaystyle|{\bm{u}}^{\mbox{\tiny{\sf
T}}}[{\bm{D}}_{n1}(\bm{\beta}_{n})-\tilde{{\bm{D}}}_{n1}(\bm{\beta}_{n})]{\bm{u}}|$
$\displaystyle\leq$
$\displaystyle||{\bm{u}}||^{2}\cdot||\widehat{{\bm{R}}}^{-1}-\tilde{{\bm{R}}}^{-1}||_{F}\cdot\lambda_{\max}({\bm{A}}_{i}(\bm{\beta}_{n}))\cdot\lambda_{\max}\Big{(}\sum_{i=1}^{n}{\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{n})\text{vec}({\bm{X}}_{i})\text{vec}^{\mbox{\tiny{\sf
T}}}({\bm{X}}_{i}){\bm{J}}(\bm{\beta}_{n})\Big{)}.$

By Condition (A4*) and (A6*), $|{\bm{u}}^{\mbox{\tiny{\sf
T}}}[{\bm{D}}_{n1}(\bm{\beta}_{n})-\tilde{{\bm{D}}}_{n1}(\bm{\beta}_{n})]{\bm{u}}|=O_{p}(\sqrt{p_{n}n})$
on the set
$\\{\bm{\beta}_{n}:||\bm{\beta}_{n}-\bm{\beta}_{0}||\leq\triangle\sqrt{p_{n}/n}\\}$.

For $i=2$, we have

$\displaystyle|{\bm{u}}^{\mbox{\tiny{\sf
T}}}[{\bm{D}}_{n2}(\bm{\beta}_{n})-\tilde{{\bm{D}}}_{n2}(\bm{\beta}_{n})]{\bm{u}}|$
$\displaystyle\leq$ $\displaystyle\frac{1}{2}|{\bm{u}}^{\mbox{\tiny{\sf
T}}}\sum_{i=1}^{n}{\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{n})\text{vec}({\bm{X}}_{i}){\bm{A}}^{1/2}_{i}(\bm{\beta}_{n})(\tilde{{\bm{R}}}^{-1}-\widehat{{\bm{R}}}^{-1}){\bm{A}}^{-3/2}_{i}(\bm{\beta}_{n}){\bf
C}_{i1}(\bm{\beta}_{n}){\bm{F}}_{i}(\bm{\beta}_{n})\text{vec}^{\mbox{\tiny{\sf
T}}}({\bm{X}}_{i}){\bm{J}}(\bm{\beta}_{n}){\bm{u}}|$
$\displaystyle+\frac{1}{2}|{\bm{u}}^{\mbox{\tiny{\sf
T}}}\sum_{i=1}^{n}{\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{n})\text{vec}({\bm{X}}_{i}){\bm{A}}^{1/2}_{i}(\bm{\beta}_{n})(\tilde{{\bm{R}}}^{-1}-\widehat{{\bm{R}}}^{-1}){\bm{A}}^{-3/2}_{i}(\bm{\beta}_{n}){\bf
C}_{i2}(\bm{\beta}_{0}){\bm{F}}_{i}(\bm{\beta}_{n})\text{vec}^{\mbox{\tiny{\sf
T}}}({\bm{X}}_{i}){\bm{J}}(\bm{\beta}_{n}){\bm{u}}|$ $\displaystyle\triangleq$
$\displaystyle J_{n1}+J_{n2},$

where

$\displaystyle{\bf C}_{i1}(\bm{\beta}_{n})$
$\displaystyle=\text{diag}\Big{(}\mu_{i1}(\bm{\beta}_{0})-\mu_{i1}(\bm{\beta}_{n}),\dots,\mu_{im}(\bm{\beta}_{0})-\mu_{im}(\bm{\beta}_{n})\Big{)},$
$\displaystyle{\bf C}_{i2}(\bm{\beta}_{0})$
$\displaystyle=\text{diag}\Big{(}Y_{i1}-\mu_{i1}(\bm{\beta}_{0}),\dots,Y_{im}-\mu_{im}(\bm{\beta}_{0})\Big{)}.$

For $J_{n1}$, by Cauchy-Schwarz inequality for matrices with Frobenius norm,

$\displaystyle J_{n1}\leq$ $\displaystyle
C||\widehat{{\bm{R}}}^{-1}-\tilde{{\bm{R}}}^{-1}||_{F}\cdot\lambda_{\max}\Big{(}\sum_{i=1}^{n}{\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{n})\text{vec}({\bm{X}}_{i})\text{vec}^{\mbox{\tiny{\sf
T}}}({\bm{X}}_{i}){\bm{J}}(\bm{\beta}_{n})\Big{)}$
$\displaystyle\max_{i,j}|\mu_{ij}^{(1)}(\widetilde{\bm{\beta}}_{n})|\cdot\mu_{ij}^{(2)}(\bm{\beta}_{n})|\cdot\frac{\max_{j}\sigma_{ij}(\bm{\beta}_{n})}{\min_{j}\sigma_{ij}(\bm{\beta}_{n})},$

where $\widetilde{\bm{\beta}}_{n}$ is between $\bm{\beta}_{n}$ and
$\bm{\beta}_{0}$. By Conditions (A3*), (A4*) and (A8*),
$\sup_{||\bm{\beta}_{n}-\bm{\beta}_{0}||\leq\triangle
n^{-1/2}}J_{n1}=O_{p}(\sqrt{p_{n}n})$.

For $J_{n2}$, we decompose ${\bm{A}}^{1/2}_{i}(\bm{\beta}_{n})$ into
${\bm{A}}^{1/2}_{i}(\bm{\beta}_{0})$ and
$[{\bm{A}}^{1/2}_{i}(\bm{\beta}_{n})-{\bm{A}}^{1/2}_{i}(\bm{\beta}_{0})]$.
That is,

$\displaystyle 2J_{n2}\leq$ $\displaystyle|{\bm{u}}^{\mbox{\tiny{\sf
T}}}\sum_{i=1}^{n}{\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{n})\text{vec}({\bm{X}}_{i}){\bm{A}}^{1/2}_{i}(\bm{\beta}_{0})(\tilde{{\bm{R}}}^{-1}-\widehat{{\bm{R}}}^{-1}){\bm{A}}^{-3/2}_{i}(\bm{\beta}_{n})$
$\displaystyle\times{\bf
C}_{i2}(\bm{\beta}_{0}){\bm{F}}_{i}(\bm{\beta}_{n})\text{vec}^{\mbox{\tiny{\sf
T}}}({\bm{X}}_{i}){\bm{J}}(\bm{\beta}_{n}){\bm{u}}|$ $\displaystyle+$
$\displaystyle|\sum_{i=1}^{n}{\bm{u}}^{\mbox{\tiny{\sf
T}}}{\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{n})\text{vec}({\bm{X}}_{i})[{\bm{A}}^{1/2}_{i}(\bm{\beta}_{n})-{\bm{A}}^{1/2}_{i}(\bm{\beta}_{0})](\tilde{{\bm{R}}}^{-1}-\widehat{{\bm{R}}}^{-1}){\bm{A}}^{-3/2}_{i}(\bm{\beta}_{n})$
$\displaystyle\times{\bf
C}_{i2}(\bm{\beta}_{0}){\bm{F}}_{i}(\bm{\beta}_{n})\text{vec}^{\mbox{\tiny{\sf
T}}}({\bm{X}}_{i}){\bm{J}}(\bm{\beta}_{n}){\bm{u}}|$ $\displaystyle\triangleq$
$\displaystyle J_{n21}+J_{n22}.$

Similarly to $J_{n1}$, it can be shown
$\sup_{||\bm{\beta}_{n}-\bm{\beta}_{0}||\leq\triangle\sqrt{p_{n}/n}}J_{n21}=O_{p}(\sqrt{np_{n}})$.

For $J_{n22}$, similar to the decomposition of
${\bm{A}}^{1/2}_{i}(\bm{\beta}_{n})$, we can further decompose those terms
involving $\bm{\beta}_{n}$ into terms that only depend on $\bm{\beta}_{0}$ and
four other terms involving $\bm{\beta}_{n}$. On the set
$\\{\bm{\beta}_{n}:||\bm{\beta}_{n}-\bm{\beta}_{0}||\leq\triangle\sqrt{p_{n}/n}\\}$,
similar to $J_{n1}$, under Conditions (A1*)-(A9*), all those terms involving
$\bm{\beta}_{n}$ can be shown to be $O_{p}(\sqrt{p_{n}n})$. To complete the
evaluation of $J_{n22}$ and hence $J_{n2}$, it suffices to show

$\displaystyle|\sum_{i=1}^{n}$ $\displaystyle{\bm{u}}^{\mbox{\tiny{\sf
T}}}{\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{0})\text{vec}({\bm{X}}_{i}){\bm{A}}^{1/2}_{i}(\bm{\beta}_{0})(\tilde{{\bm{R}}}^{-1}-\widehat{{\bm{R}}}^{-1}){\bm{A}}^{-3/2}_{i}(\bm{\beta}_{0})$
$\displaystyle\times{\bf
C}_{i2}(\bm{\beta}_{0}){\bm{F}}_{i}(\bm{\beta}_{0})\text{vec}^{\mbox{\tiny{\sf
T}}}({\bm{X}}_{i}){\bm{J}}(\bm{\beta}_{0}){\bm{u}}|=O_{p}(\sqrt{p_{n}n}).$
(10)

Denote $L_{n}(\bm{\beta}_{0})$ the left side of (Proof of Lemma 3.). Recall
that
$\epsilon_{ij}(\bm{\beta}_{0})=\sigma_{ij}^{-1}(\bm{\beta}_{0})(Y_{ij}-\mu_{ij}(\bm{\beta}_{0}))$.
We have

$\displaystyle\mathbb{E}[||L_{n}(\bm{\beta}_{0})||^{2}]=\text{Tr}[\mathbb{E}(L_{n}(\bm{\beta}_{0})^{\mbox{\tiny{\sf
T}}}L_{n}(\bm{\beta}_{0}))]$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{n}\sum_{j=1}^{m}\sum_{k=1}^{m}\mathbb{E}[\epsilon_{ij}\epsilon_{ik}]\text{Tr}\Big{[}{\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{0})\text{vec}({\bm{X}}_{i}){\bm{A}}^{1/2}_{i}(\bm{\beta}_{0})(\tilde{{\bm{R}}}^{-1}-\widehat{{\bm{R}}}^{-1}){\bm{A}}^{-3/2}_{i}(\bm{\beta}_{0}){\bf
e}_{j}{\bf e}_{j}^{\mbox{\tiny{\sf T}}}{\bm{F}}_{i}(\bm{\beta}_{0})$
$\displaystyle\cdot\text{vec}^{\mbox{\tiny{\sf
T}}}({\bm{X}}_{i}){\bm{J}}(\bm{\beta}_{0}){\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{0})\text{vec}({\bm{X}}_{i}){\bf e}_{k}{\bf
e}_{k}^{\mbox{\tiny{\sf
T}}}{\bm{A}}^{-3/2}_{i}(\bm{\beta}_{0})(\tilde{{\bm{R}}}^{-1}-\widehat{{\bm{R}}}^{-1}){\bm{A}}^{1/2}_{i}(\bm{\beta}_{0})\text{vec}^{\mbox{\tiny{\sf
T}}}({\bm{X}}_{i}){\bm{J}}(\bm{\beta}_{0})$ $\displaystyle\leq$ $\displaystyle
C\sum_{i=1}^{n}\sum_{j=1}^{m}\sum_{k=1}^{m}||{\bf e}_{j}^{\mbox{\tiny{\sf
T}}}{\bm{F}}_{i}(\bm{\beta}_{0})\text{vec}^{\mbox{\tiny{\sf
T}}}({\bm{X}}_{i}){\bm{J}}(\bm{\beta}_{0})||\cdot||{\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{0})\text{vec}({\bm{X}}_{i}){\bm{F}}_{i}(\bm{\beta}_{0}){\bf
e}_{k}||$ $\displaystyle\cdot||{\bf e}_{k}^{\mbox{\tiny{\sf
T}}}{\bm{A}}^{-3/2}_{i}(\bm{\beta}_{0})(\tilde{{\bm{R}}}^{-1}-\widehat{{\bm{R}}}^{-1}){\bm{A}}^{1/2}_{i}(\bm{\beta}_{0})\text{vec}^{\mbox{\tiny{\sf
T}}}({\bm{X}}_{i}){\bm{J}}(\bm{\beta}_{0})||$
$\displaystyle\cdot||{\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{0})\text{vec}({\bm{X}}_{i}){\bm{A}}^{1/2}_{i}(\bm{\beta}_{0})(\tilde{{\bm{R}}}^{-1}-\widehat{{\bm{R}}}^{-1}){\bm{A}}^{-3/2}_{i}(\bm{\beta}_{0}){\bf
e}_{j}||.$

By Conditions (A1*), (A2*) and (A4*)-(A7*),
$\mathbb{E}[||L_{n}(\bm{\beta}_{0})||^{2}]=O(np_{n})$. This implies
$J_{n22}=O_{p}(\sqrt{p_{n}n})$ and hence

$\sup_{||\bm{\beta}_{n}-\bm{\beta}_{0}||\leq\triangle\sqrt{p_{n}/n}}\sup_{{\bm{u}}}|{\bm{u}}^{\mbox{\tiny{\sf
T}}}[{\bm{D}}_{n2}(\bm{\beta}_{n})-\tilde{{\bm{D}}}_{n2}(\bm{\beta}_{n})]{\bm{u}}|=O_{p}(\sqrt{p_{n}n}).$

Using similar decompositions, we can verify the results for ${\bm{D}}_{n3}$
and ${\bm{D}}_{n4}$, which completes the proof. ∎

Based on Lemma 3, we can further approximate
${\bf\tilde{D}}_{n}(\bm{\beta}_{n})$ by ${\bf\tilde{D}}_{n1}(\bm{\beta}_{n})$,
which is easier to evaluate. Lemma 4 provides this approximation.

######  Lemma 4.

Under Conditions (A1*)-(A9*), for any $\triangle>0$ and
${\bm{u}}\in\mathrm{I\\!R}\mathit{{}^{R\sum_{d=1}^{D}p_{d}}}$ such that
$||{\bm{u}}||=1$,

$\displaystyle\sup_{||\bm{\beta}_{n}-\bm{\beta}_{0}||=\triangle\sqrt{p_{n}/n}}\sup_{{\bm{u}}}|{\bm{u}}^{\mbox{\tiny{\sf
T}}}[\tilde{{\bm{D}}}_{n}(\bm{\beta}_{n})-\tilde{{\bm{D}}}_{n1}(\bm{\beta}_{n})]{\bm{u}}|=O_{p}(n^{1/2}p_{n}),$
(11)
$\displaystyle\sup_{||\bm{\beta}_{n}-\bm{\beta}_{0}||=\triangle\sqrt{p_{n}/n}}\sup_{{\bm{u}}}|{\bm{u}}^{\mbox{\tiny{\sf
T}}}[\tilde{{\bm{D}}}_{n1}(\bm{\beta}_{0})-\tilde{{\bm{D}}}_{n1}(\bm{\beta}_{n})]{\bm{u}}|=O_{p}(n^{1/2}p_{n}).$
(12)

###### Proof of Lemma 4.

To prove (11), it is sufficient to show, for $i=2,3,4$,

$\sup_{||\bm{\beta}_{n}-\bm{\beta}_{0}||=\triangle\sqrt{p_{n}/n}}\sup_{{\bm{u}}}|{\bm{u}}^{\mbox{\tiny{\sf
T}}}\tilde{{\bm{D}}}_{ni}(\bm{\beta}_{n}){\bm{u}}|=O_{p}(n^{1/2}p_{n}).$

For $\tilde{{\bm{D}}}_{n2}(\bm{\beta}_{n})$, we have the decomposition

$\tilde{{\bm{D}}}_{n2}(\bm{\beta}_{n})=\tilde{{\bm{D}}}_{n2}(\bm{\beta}_{0})+\sum_{k=1}^{6}J_{n6}(\bm{\beta}_{n}),\\\
$

where

$\displaystyle J_{n1}$
$\displaystyle=\frac{1}{2}\sum_{i=1}^{n}[{\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{n})-{\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{0})]\text{vec}({\bm{X}}_{i}){\bm{A}}^{1/2}_{i}(\bm{\beta}_{0})\tilde{{\bm{R}}}^{-1}{\bm{A}}^{-3/2}_{i}(\bm{\beta}_{0}){\bf
C}_{i}(\bm{\beta}_{0}){\bm{F}}_{i}(\bm{\beta}_{0})\text{vec}^{\mbox{\tiny{\sf
T}}}({\bm{X}}_{i}){\bm{J}}(\bm{\beta}_{0}),$ $\displaystyle J_{n2}$
$\displaystyle=\frac{1}{2}\sum_{i=1}^{n}{\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{n})\text{vec}({\bm{X}}_{i})[{\bm{A}}^{1/2}_{i}(\bm{\beta}_{n})-{\bm{A}}^{1/2}_{i}(\bm{\beta}_{0})]\tilde{{\bm{R}}}^{-1}{\bm{A}}^{-3/2}_{i}(\bm{\beta}_{0}){\bf
C}_{i}(\bm{\beta}_{0}){\bm{F}}_{i}(\bm{\beta}_{0})\text{vec}^{\mbox{\tiny{\sf
T}}}({\bm{X}}_{i}){\bm{J}}(\bm{\beta}_{0}),$ $\displaystyle J_{n3}$
$\displaystyle=\frac{1}{2}\sum_{i=1}^{n}{\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{n})\text{vec}({\bm{X}}_{i}){\bm{A}}^{1/2}_{i}(\bm{\beta}_{n})\tilde{{\bm{R}}}^{-1}[{\bm{A}}^{-3/2}_{i}(\bm{\beta}_{n})-{\bm{A}}^{-3/2}_{i}(\bm{\beta}_{0})]{\bf
C}_{i}(\bm{\beta}_{0}){\bm{F}}_{i}(\bm{\beta}_{0})\text{vec}^{\mbox{\tiny{\sf
T}}}({\bm{X}}_{i}){\bm{J}}(\bm{\beta}_{0}),$ $\displaystyle J_{n4}$
$\displaystyle=\frac{1}{2}\sum_{i=1}^{n}{\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{n})\text{vec}({\bm{X}}_{i}){\bm{A}}^{1/2}_{i}(\bm{\beta}_{n})\tilde{{\bm{R}}}^{-1}{\bm{A}}^{-3/2}_{i}(\bm{\beta}_{0})[{\bf
C}_{i}(\bm{\beta}_{n})-{\bf
C}_{i}(\bm{\beta}_{0})]{\bm{F}}_{i}(\bm{\beta}_{0})\text{vec}^{\mbox{\tiny{\sf
T}}}({\bm{X}}_{i}){\bm{J}}(\bm{\beta}_{0}),$ $\displaystyle J_{n5}$
$\displaystyle=\frac{1}{2}\sum_{i=1}^{n}{\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{n})\text{vec}({\bm{X}}_{i}){\bm{A}}^{1/2}_{i}(\bm{\beta}_{n})\tilde{{\bm{R}}}^{-1}{\bm{A}}^{-3/2}_{i}(\bm{\beta}_{n}){\bf
C}_{i}(\bm{\beta}_{n})[{\bm{F}}_{i}(\bm{\beta}_{n})-{\bm{F}}_{i}(\bm{\beta}_{0})]\text{vec}^{\mbox{\tiny{\sf
T}}}({\bm{X}}_{i}){\bm{J}}(\bm{\beta}_{0}),$ $\displaystyle J_{n6}$
$\displaystyle=\frac{1}{2}\sum_{i=1}^{n}{\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{n})\text{vec}({\bm{X}}_{i}){\bm{A}}^{1/2}_{i}(\bm{\beta}_{n})\tilde{{\bm{R}}}^{-1}{\bm{A}}^{-3/2}_{i}(\bm{\beta}_{n}){\bf
C}_{i}(\bm{\beta}_{n}){\bm{F}}_{i}(\bm{\beta}_{n})\text{vec}^{\mbox{\tiny{\sf
T}}}({\bm{X}}_{i})[{\bm{J}}(\bm{\beta}_{n})-{\bm{J}}(\bm{\beta}_{0})].$

Using the same techniques as in Lemma 3, it can be shown that
$\tilde{{\bm{D}}}_{n2}(\bm{\beta}_{0})$ and $J_{ni}$ are
$O_{p}(n^{1/2}p_{n})$, $i=1,\dots,6$. We can prove
$\tilde{{\bm{D}}}_{n3}(\bm{\beta}_{0})$ and
$\tilde{{\bm{D}}}_{n4}(\bm{\beta}_{0})$ in the same way, which completes the
proof of (11).

To prove (12), note that

$\displaystyle|{\bm{u}}^{\mbox{\tiny{\sf
T}}}[\tilde{{\bm{D}}}_{n1}(\bm{\beta}_{0})-\tilde{{\bm{D}}}_{n1}(\bm{\beta}_{n})]{\bm{u}}|$
$\displaystyle\leq$ $\displaystyle|{\bm{u}}^{\mbox{\tiny{\sf
T}}}[{\bm{J}}^{\mbox{\tiny{\sf T}}}(\bm{\beta}_{0})-{\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{n})]\text{vec}({\bm{X}}_{i}){\bm{A}}^{1/2}_{i}(\bm{\beta}_{n})\tilde{{\bm{R}}}^{-1}{\bm{A}}^{1/2}_{i}(\bm{\beta}_{n})\text{vec}^{\mbox{\tiny{\sf
T}}}({\bm{X}}_{i}){\bm{J}}(\bm{\beta}_{n}){\bm{u}}|$
$\displaystyle+|{\bm{u}}^{\mbox{\tiny{\sf T}}}{\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{0})\text{vec}({\bm{X}}_{i})[{\bm{A}}^{1/2}_{i}(\bm{\beta}_{0})-{\bm{A}}^{1/2}_{i}(\bm{\beta}_{n})]\tilde{{\bm{R}}}^{-1}{\bm{A}}^{1/2}_{i}(\bm{\beta}_{n})\text{vec}^{\mbox{\tiny{\sf
T}}}({\bm{X}}_{i}){\bm{J}}(\bm{\beta}_{n}){\bm{u}}|$
$\displaystyle+|{\bm{u}}^{\mbox{\tiny{\sf T}}}{\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{0})\text{vec}({\bm{X}}_{i}){\bm{A}}^{1/2}_{i}(\bm{\beta}_{0})\tilde{{\bm{R}}}^{-1}[{\bm{A}}^{1/2}_{i}(\bm{\beta}_{0})-{\bm{A}}^{1/2}_{i}(\bm{\beta}_{n})]\text{vec}^{\mbox{\tiny{\sf
T}}}({\bm{X}}_{i}){\bm{J}}(\bm{\beta}_{n}){\bm{u}}|$
$\displaystyle+|{\bm{u}}^{\mbox{\tiny{\sf T}}}{\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{0})\text{vec}({\bm{X}}_{i}){\bm{A}}^{1/2}_{i}(\bm{\beta}_{0})\tilde{{\bm{R}}}^{-1}{\bm{A}}^{1/2}_{i}(\bm{\beta}_{0})\text{vec}^{\mbox{\tiny{\sf
T}}}({\bm{X}}_{i})[{\bm{J}}(\bm{\beta}_{0})-{\bm{J}}(\bm{\beta}_{n})]{\bm{u}}|.$

The rest of the proof is similar to Lemma 3 and thus is omitted here. ∎

Proof of Theorem 3

###### Proof.

Wang, (2011) gave a sufficient condition for the existence and consistency of
a sequence of root $\widehat{\bm{\beta}}_{n}$of
${\bm{s}}_{n}(\bm{\beta}_{n})=0$, namely,

$P(\sup_{||\bm{\beta}_{n}-\bm{\beta}_{0}||=\triangle\sqrt{p_{n}/n}}(\bm{\beta}_{n}-\bm{\beta}_{0})^{\mbox{\tiny{\sf
T}}}{\bm{s}}_{n}(\bm{\beta}_{n})<0)\geq 1-\epsilon$ (13)

with $\forall\epsilon>0$ and a constant $\triangle>0$. To verify (13), the
main idea is to approximate ${\bm{s}}_{n}(\bm{\beta}_{n})$ by
$\tilde{{\bm{s}}}_{n}(\bm{\beta}_{n})$, whose moments are easier to evaluate.

By direct calculation,

$\displaystyle(\bm{\beta}_{n}-\bm{\beta}_{0})^{\mbox{\tiny{\sf
T}}}{\bm{s}}_{n}(\bm{\beta}_{n})$ $\displaystyle=$
$\displaystyle(\bm{\beta}_{n}-\bm{\beta}_{0})^{\mbox{\tiny{\sf
T}}}{\bm{s}}_{n}(\bm{\beta}_{0})-(\bm{\beta}_{n}-\bm{\beta}_{0})^{\mbox{\tiny{\sf
T}}}{\bm{D}}_{n}(\bm{\beta}^{*}_{n})(\bm{\beta}_{n}-\bm{\beta}_{0})$
$\displaystyle\triangleq$ $\displaystyle I_{n1}+I_{n2},$

where $\bm{\beta}^{*}_{n}=t\bm{\beta}_{n}+(1-t)\bm{\beta}^{0}$ for some
$0<t<1$. Further decompose $I_{n1}$ into

$\displaystyle I_{n1}$
$\displaystyle=(\bm{\beta}_{n}-\bm{\beta}_{0})^{\mbox{\tiny{\sf
T}}}\tilde{{\bm{s}}}_{n}(\bm{\beta}_{0})+(\bm{\beta}_{n}-\bm{\beta}_{0})^{\mbox{\tiny{\sf
T}}}[{\bm{s}}_{n}(\bm{\beta}_{0})-\tilde{{\bm{s}}}_{n}(\bm{\beta}_{0})]$
$\displaystyle\triangleq I_{n11}+I_{n12}.$

Note that
$I_{n11}\leq\triangle\sqrt{p_{n}/n}\cdot||\tilde{{\bm{s}}}_{n}(\bm{\beta}_{0})||$,
where

$\displaystyle\mathbb{E}[||\tilde{{\bm{s}}}_{n}(\bm{\beta}_{0})||^{2}]$
$\displaystyle=$
$\displaystyle\mathbb{E}\Big{\\{}\sum_{i=1}^{n}\bm{\epsilon}_{i}^{\mbox{\tiny{\sf
T}}}\tilde{{\bm{R}}}^{-1}{\bm{A}}^{1/2}_{i}(\bm{\beta}_{0})\text{vec}^{\mbox{\tiny{\sf
T}}}({\bm{X}}_{i}){\bm{J}}(\bm{\beta}_{n}){\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{n})\text{vec}({\bm{X}}_{i}){\bm{A}}^{1/2}_{i}(\bm{\beta}_{0})\tilde{{\bm{R}}}^{-1}\bm{\epsilon}_{i}\Big{\\}}$
$\displaystyle\leq$ $\displaystyle
C\cdot\text{Tr}\Big{(}{\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{n})\text{vec}({\bm{X}}_{i})\text{vec}^{\mbox{\tiny{\sf
T}}}({\bm{X}}_{i}){\bm{J}}(\bm{\beta}_{n})\Big{)}$ $\displaystyle=$
$\displaystyle C\sum_{j=1}^{m}\cdot\text{Tr}\Big{(}{\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{n})\text{vec}({\bm{X}}_{ij})\text{vec}^{\mbox{\tiny{\sf
T}}}({\bm{X}}_{ij}){\bm{J}}(\bm{\beta}_{n})\Big{)}=O(np_{n})$

for some constant $C>0$. This implies that
$I_{n11}=\triangle\sqrt{p_{n}/n}O_{p}(\sqrt{np_{n}})=\triangle O_{p}(p_{n})$.
For $I_{n12}$, by Lemma 2,

$I_{n12}\leq||\bm{\beta}_{n}-\bm{\beta}_{0}||\cdot||{\bm{s}}_{n}(\bm{\beta}_{0})-\tilde{{\bm{s}}}_{n}(\bm{\beta}_{0})||=o_{p}(p_{n}).$

Therefore, $I_{n1}$ is dominated in probability by $I_{n11}$.

For $I_{n2}$ ,we decompose it into

$\displaystyle I_{n2}=$
$\displaystyle-(\bm{\beta}_{n}-\bm{\beta}_{0})^{\mbox{\tiny{\sf
T}}}\tilde{{\bm{D}}}_{n}(\bm{\beta}^{*}_{n})(\bm{\beta}_{n}-\bm{\beta}_{0})$
$\displaystyle-(\bm{\beta}_{n}-\bm{\beta}_{0})^{\mbox{\tiny{\sf
T}}}[{\bm{D}}_{n}(\bm{\beta}^{*}_{n})-\tilde{{\bm{D}}}_{n}(\bm{\beta}^{*}_{n})](\bm{\beta}_{n}-\bm{\beta}_{0})$
$\displaystyle\triangleq$ $\displaystyle I_{n21}+I_{n22}.$

By Lemma 3, it can be easily checked that $I_{n22}=o_{p}(p_{n})$. Next for
$I_{n21}$,

$\displaystyle I_{n21}=$
$\displaystyle-(\bm{\beta}_{n}-\bm{\beta}_{0})^{\mbox{\tiny{\sf
T}}}\tilde{{\bm{D}}}_{n1}(\bm{\beta}_{0})(\bm{\beta}_{n}-\bm{\beta}_{0})$
$\displaystyle-(\bm{\beta}_{n}-\bm{\beta}_{0})^{\mbox{\tiny{\sf
T}}}[\tilde{{\bm{D}}}_{n1}(\bm{\beta}^{*}_{n})-\tilde{{\bm{D}}}_{n1}(\bm{\beta}_{0})](\bm{\beta}_{n}-\bm{\beta}_{0})$
$\displaystyle-(\bm{\beta}_{n}-\bm{\beta}_{0})^{\mbox{\tiny{\sf
T}}}[\tilde{{\bm{D}}}_{n}(\bm{\beta}^{*}_{n})-\tilde{{\bm{D}}}_{n1}(\bm{\beta}^{*}_{n})](\bm{\beta}_{n}-\bm{\beta}_{0})$
$\displaystyle\triangleq$ $\displaystyle I^{1}_{n21}+I^{2}_{n21}+I^{3}_{n21}.$

We next show that $I_{n21}$ is dominated in probability by $I^{1}_{n21}$. Note
that by Condition (A3*), (A4*) and (A8*),

$\displaystyle I^{1}_{n21}=$
$\displaystyle-(\bm{\beta}_{n}-\bm{\beta}_{0})^{\mbox{\tiny{\sf
T}}}\Big{[}\sum_{i=1}^{n}{\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{n})\text{vec}({\bm{X}}_{i}){\bm{A}}^{1/2}_{i}(\bm{\beta}_{n})\tilde{{\bm{R}}}^{-1}{\bm{A}}^{1/2}_{i}(\bm{\beta}_{n})\text{vec}^{\mbox{\tiny{\sf
T}}}({\bm{X}}_{i}){\bm{J}}(\bm{\beta}_{n})\Big{]}(\bm{\beta}_{n}-\bm{\beta}_{0})$
$\displaystyle\leq$
$\displaystyle-n^{-1}\triangle^{2}\min_{i}\lambda_{\min}({\bm{A}}_{i}(\bm{\beta}_{n}))\lambda_{\min}\Big{(}\sum_{i=1}^{n}{\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{n})\text{vec}({\bm{X}}_{i})\text{vec}^{\mbox{\tiny{\sf
T}}}({\bm{X}}_{i}){\bm{J}}(\bm{\beta}_{n})\Big{)}\lambda_{\min}(\tilde{{\bm{R}}}^{-1})$
$\displaystyle\leq$ $\displaystyle-C\triangle^{2}p_{n},$

for some constant $C>0$. By Lemma 4, it can be checked directly that both
$I^{2}_{n21}$ and $I^{3}_{n21}$ are $o_{p}(p_{n})$.

Therefore, the sign of $(\bm{\beta}_{n}-\bm{\beta}_{0})^{\mbox{\tiny{\sf
T}}}{\bm{s}}_{n}(\bm{\beta}_{n})$ is determined by in probability by
$I_{n11}+I^{1}_{n21}$ and is negative for sufficiently large $\triangle$,
which completes the proof.∎

Proof of Theorem 4

###### Proof.

We first show that the normalized $\tilde{{\bm{s}}}_{n}(\bm{\beta}_{0})$ has
an asymptotic normal distribution. That is, for any
${\bm{b}}\in\mathrm{I\\!R}\mathit{{}^{R\sum_{d=1}^{D}p_{d}}}$ such that
$||{\bm{b}}||=1$,

${\bm{b}}^{\mbox{\tiny{\sf
T}}}\tilde{{\bm{M}}}^{-1/2}_{n}(\bm{\beta}_{0})\tilde{{\bm{s}}}_{n}(\bm{\beta}_{0})\to
N(0,1),$ (14)

where
$\tilde{{\bm{M}}}_{n}(\bm{\beta}_{0})=\text{Var}(\tilde{{\bm{s}}}_{n}(\bm{\beta}_{0}))$.

Denote ${\bm{b}}^{\mbox{\tiny{\sf
T}}}\tilde{{\bm{M}}}^{-1/2}_{n}(\bm{\beta}_{0})\tilde{{\bm{s}}}_{n}(\bm{\beta}_{0})=\sum_{i=1}^{n}Z_{ni}$,
where

$Z_{ni}={\bm{b}}^{\mbox{\tiny{\sf
T}}}\tilde{{\bm{M}}}^{-1/2}_{n}(\bm{\beta}_{0}){\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{0})\text{vec}({\bm{X}}_{i}){\bm{A}}^{1/2}_{i}(\bm{\beta}_{0})\tilde{{\bm{R}}}^{-1}\bm{\epsilon}_{i}(\bm{\beta}_{0}),$

where
$\bm{\epsilon}_{i}(\bm{\beta}_{0})={\bm{A}}^{-1/2}_{i}(\bm{\beta}_{0})({\bf
Y}_{i}-\bm{\mu}_{i}(\bm{\beta}_{0}))$. Note that $\mathbb{E}(Z_{ni})=0$,
$\text{Var}(\sum_{i=1}^{n}Z_{ni})=1$. To prove (14), it suffices to check the
Lyapunov condition. That is, for some $\delta>0$,

$\sum_{i=1}^{n}\mathbb{E}\Big{(}|Z_{ni}|^{2+\delta}\Big{)}\to 0,$

as $n\to\infty$. By Cauchy-Schwarz inequality,

$\displaystyle
Z_{ni}^{2}\leq\lambda_{\max}(\tilde{{\bm{R}}}^{-2})\lambda_{\max}({\bm{A}}_{i}(\bm{\beta}_{0}))||\bm{\epsilon}_{i}(\bm{\beta}_{0})||^{2}\gamma_{ni},$

where $\gamma_{ni}\triangleq{\bm{b}}^{\mbox{\tiny{\sf
T}}}\tilde{{\bm{M}}}^{-1/2}_{n}(\bm{\beta}_{0}){\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{0})\text{vec}({\bm{X}}_{i})\text{vec}^{\mbox{\tiny{\sf
T}}}({\bm{X}}_{i}){\bm{J}}(\bm{\beta}_{0})\tilde{{\bm{M}}}^{-1/2}_{n}(\bm{\beta}_{0}){\bm{b}}$.
To evaluate $\max_{1\leq i\leq n}\gamma_{ni}$, we need to evaluate
$\lambda^{-1}_{\min}(\tilde{{\bm{M}}}_{n}(\bm{\beta}_{0}))$. Note that

$\displaystyle{\bm{b}}^{\mbox{\tiny{\sf
T}}}\tilde{{\bm{M}}}_{n}(\bm{\beta}_{0}){\bm{b}}\geq$ $\displaystyle
C{\bm{b}}^{\mbox{\tiny{\sf T}}}\Big{(}\sum_{i=1}^{n}{\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{0})\text{vec}({\bm{X}}_{i})\text{vec}^{\mbox{\tiny{\sf
T}}}({\bm{X}}_{i}){\bm{J}}(\bm{\beta}_{0})\Big{)}{\bm{b}}$ $\displaystyle\geq$
$\displaystyle C\lambda_{\min}\Big{(}\sum_{i=1}^{n}{\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{0})\text{vec}({\bm{X}}_{i})\text{vec}^{\mbox{\tiny{\sf
T}}}({\bm{X}}_{i}){\bm{J}}(\bm{\beta}_{0})\Big{)}.$

By Condition (A1*) and (A3*),
$\lambda^{-1}_{\min}(\tilde{{\bm{M}}}_{n}(\bm{\beta}_{0}))=O(n^{-1})$ and
hence $\max_{1\leq i\leq n}\gamma_{ni}=o(1)$.

It follows that for any $\delta>0$,

$\displaystyle\sum_{i=1}^{n}\mathbb{E}\Big{(}|Z_{ni}|^{2+\delta}\Big{)}$
$\displaystyle\leq\sum_{i=1}^{n}\mathbb{E}\Big{(}C^{1+\delta/2}\gamma_{ni}^{1+\delta/2}||\bm{\epsilon}_{i}(\bm{\beta}_{0})||^{2+\delta}\Big{)}$
$\displaystyle\leq C(\max_{1\leq i\leq
n}\gamma_{ni})^{\delta/2}\sum_{i=1}^{n}{\bm{b}}^{\mbox{\tiny{\sf
T}}}\tilde{{\bm{M}}}^{-1/2}_{n}(\bm{\beta}_{0}){\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{0})\text{vec}({\bm{X}}_{i})\text{vec}^{\mbox{\tiny{\sf
T}}}({\bm{X}}_{i}){\bm{J}}(\bm{\beta}_{0})\tilde{{\bm{M}}}^{-1/2}_{n}(\bm{\beta}_{0}){\bm{b}}$
$\displaystyle\leq C(\max_{1\leq i\leq
n}\gamma_{ni})^{\delta/2}\lambda_{\max}\big{(}\sum_{i=1}^{n}{\bm{J}}^{\mbox{\tiny{\sf
T}}}(\bm{\beta}_{0})\text{vec}({\bm{X}}_{i})\text{vec}^{\mbox{\tiny{\sf
T}}}({\bm{X}}_{i}){\bm{J}}(\bm{\beta}_{0})\big{)}\lambda^{-1}_{\min}\big{(}\tilde{{\bm{M}}}_{n}(\bm{\beta}_{0})\big{)}$
$\displaystyle=o(1)O(n)O(n^{-1})=o(1),$

which completes the proof of (14).

To prove Theorem 2, note that by the fact
${\bm{s}}_{n}(\widehat{\bm{\beta}}_{n})=0$, we have
${\bm{s}}_{n}(\bm{\beta}_{0})={\bm{D}}_{n}(\bm{\beta}^{*}_{n})(\widehat{\bm{\beta}}_{n}-\bm{\beta}_{0})$
for some $\bm{\beta}^{*}_{n}$ between $\widehat{\bm{\beta}}_{n}$ and
$\bm{\beta}_{0}$. Hence,

$\displaystyle{\bm{b}}^{\mbox{\tiny{\sf
T}}}\tilde{{\bm{M}}}^{-1/2}_{n}(\bm{\beta}_{0})\tilde{{\bm{s}}}_{n}(\bm{\beta}_{0})$
$\displaystyle=$ $\displaystyle{\bm{b}}^{\mbox{\tiny{\sf
T}}}\tilde{{\bm{M}}}^{-1/2}_{n}(\bm{\beta}_{0})\tilde{{\bm{D}}}_{n1}(\bm{\beta}_{0})(\widehat{\bm{\beta}}_{n}-\bm{\beta}_{0})$
$\displaystyle+{\bm{b}}^{\mbox{\tiny{\sf
T}}}\tilde{{\bm{M}}}^{-1/2}_{n}(\bm{\beta}_{0})[{\bm{D}}_{n}(\bm{\beta}^{*}_{n})-\tilde{{\bm{D}}}_{n1}(\bm{\beta}_{0})](\widehat{\bm{\beta}}_{n}-\bm{\beta}_{0})$
$\displaystyle+{\bm{b}}^{\mbox{\tiny{\sf
T}}}\tilde{{\bm{M}}}^{-1/2}_{n}(\bm{\beta}_{0})[\tilde{{\bm{s}}}_{n}(\bm{\beta}_{0})-{\bm{s}}_{n}(\bm{\beta}_{0})]$
$\displaystyle=$ $\displaystyle
J_{n1}+J_{n2}(\bm{\beta}^{*}_{n})+J_{n3}(\bm{\beta}_{0}).$

By (14), it is sufficient to prove that both
$\sup_{||\bm{\beta}_{n}-\bm{\beta}_{0}||\leq\triangle\sqrt{p_{n}/n}}|J_{n2}(\bm{\beta}_{n})|$
and $|J_{n3}(\bm{\beta}_{0})|$ are $o_{p}(1)$.

For $J_{n3}$, recall that
$||\tilde{{\bm{s}}}_{n}(\bm{\beta}_{0})-{\bm{s}}_{n}(\bm{\beta}_{0})||=O_{p}(1)$
from Lemma 2. Using the previous result that
$\lambda^{-1}_{\min}(\tilde{{\bm{M}}}_{n}(\bm{\beta}_{0}))=O(n^{-1})$, it can
be easily checked that $J_{n3}^{2}=o_{p}(1)$ and hence $|J_{n3}|=o_{p}(1)$.

For $J_{n2}$, we have

$\displaystyle\sup_{||\bm{\beta}_{n}-\bm{\beta}_{0}||\leq\triangle\sqrt{p_{n}/n}}|J_{n2}(\bm{\beta}_{n})|$
$\displaystyle\leq$
$\displaystyle\sup_{||\bm{\beta}_{n}-\bm{\beta}_{0}||\leq\triangle\sqrt{p_{n}/n}}{\bm{b}}^{\mbox{\tiny{\sf
T}}}\tilde{{\bm{M}}}^{-1/2}_{n}(\bm{\beta}_{0})[{\bm{D}}_{n}(\bm{\beta}_{n})-\tilde{{\bm{D}}}_{n}(\bm{\beta}_{n})](\widehat{\bm{\beta}}_{n}-\bm{\beta}_{0})$
$\displaystyle+\sup_{||\bm{\beta}_{n}-\bm{\beta}_{0}||\leq\triangle\sqrt{p_{n}/n}}{\bm{b}}^{\mbox{\tiny{\sf
T}}}\tilde{{\bm{M}}}^{-1/2}_{n}(\bm{\beta}_{0})[\tilde{{\bm{D}}}_{n}(\bm{\beta}_{n})-\tilde{{\bm{D}}}_{n1}(\bm{\beta}_{n})](\widehat{\bm{\beta}}_{n}-\bm{\beta}_{0})$
$\displaystyle+\sup_{||\bm{\beta}_{n}-\bm{\beta}_{0}||\leq\triangle\sqrt{p_{n}/n}}{\bm{b}}^{\mbox{\tiny{\sf
T}}}\tilde{{\bm{M}}}^{-1/2}_{n}(\bm{\beta}_{0})[\tilde{{\bm{D}}}_{n1}(\bm{\beta}_{n})-\tilde{{\bm{D}}}_{n1}(\bm{\beta}_{0})](\widehat{\bm{\beta}}_{n}-\bm{\beta}_{0})$
$\displaystyle\triangleq$ $\displaystyle I_{n1}+I_{n2}+I_{n3}.$

By Theorem 1, Lemma 4 and the fact
$\lambda^{-1}_{\min}(\tilde{{\bm{M}}}_{n}(\bm{\beta}_{0}))=O(n^{-1})$, it can
be seen that

$\displaystyle I_{n1}\leq$ $\displaystyle
C\cdot\lambda_{\max}({\bm{D}}_{n}(\bm{\beta}_{n})-\tilde{{\bm{D}}}_{n}(\bm{\beta}_{n}))\lambda^{-1/2}_{\min}(\tilde{{\bm{M}}}_{n}(\bm{\beta}_{0}))\sqrt{p_{n}/n}$
$\displaystyle=$ $\displaystyle
O_{p}(\sqrt{n}p_{n})O(n^{-1/2})O_{p}(\sqrt{p_{n}/n})=O_{p}(p_{n}^{3/2}n^{-1/2}).$

Using the stronger assumption that $p=o(n^{-1/3})$, $I_{n1}=o_{p}(1)$.
Similarly, by Lemma 4, we have $I_{n2}=o_{p}(1)$ and $I_{n3}=o_{p}(1)$.
Therefore $J_{n1}$ has the same asymptotic distribution as in (14), which
completes the proof. ∎