# Heterogeneous Overdispersed Count Data Regressions via Double Penalized
Estimations

Shaomin Li111Center for Statistics and Data Science, Beijing Normal University
at Zhuhai. Email<EMAIL_ADDRESS>(Shaomin Li). This work is supported in
part by National Natural Science Foundation of China Grants 12101056., Haoyu
Wei222Corresponding Author. Guanghua School of Management, Peking University.
Email<EMAIL_ADDRESS>(Haoyu Wei). The first two authors are co-first
authors, in alphabetical order., Xiaoyu Lei333Department of Statistics,
University of Chicago. Email<EMAIL_ADDRESS>(Xiaoyu Lei)

###### Abstract

This paper studies the non-asymptotic merits of the double $\ell_{1}$-penalty
for heterogeneous overdispersed count data via negative binomial regressions.
Under the restricted eigenvalue conditions, we prove the oracle inequalities
for Lasso estimators of two partial regression coefficients for the first
time, using concentration inequalities of empirical processes. Furthermore,
derived from the oracle inequalities, the consistency and convergence rate for
the estimators are the theoretical guarantees for further statistical
inference. Finally, both simulations and a real data analysis demonstrate that
the new methods are effective.

Keywords: Negative binomial regressions; Heterogeneous count data regression;
Estimation of dispersion parameter; Oracle inequalities.

##  1 Introduction

With the advance of modern data collection techniques, scientists and
engineers generate or get access to a massive number of variables in their
experiments, and challenges to traditional statistical methods and theories
have come. One crucial aspect is the high-dimensional setting, where the
number of covariates can be comparable to or greater than the sample size. In
high-dimensional settings, the asymptotical results for the estimator are
intractable. For example, the maximum likelihood estimator (MLE) in classical
multivariate statistics is easy to obtain since the Hessian matrices in
optimizations are invertible, and classical Newton algorithms also perform
well. When the dimension is high, however, the optimizations for the MLE
method lead to many un-meaningful solutions, so the specific constraint for
variables in the optimizations is essential. The Lasso-regularization method
is one of the constrained least-squares methods widely applied in the high-
dimensional parameter estimation problem; see Tibshirani (1996).

In many scientific fields such as biomedical science, ecology, and economics,
experimental and observational studies often yield count data, a type of data
in which the observations can take only the non-negative integer values. The
Poisson regression models are commonly used for count data. However, it needs
a restrictive assumption that the variance equals the mean. For many count
data, the variance is often larger than the mean (Dai et al., 2013; Zhang et
al., 2018), which is called over-dispersion. Since the Poisson regression
model is invalid under the over-dispersion case, a more general and flexible
regression model, the negative binomial regression (NBR), has attracted lots
of research attention and become popular in analyzing count data (Xie and
Xiao, 2020). Recently, there has been much research on the high-dimensional
NBR model, such as Qiu et al. (2018); Weißbach and Radloff (2020); Tian et al.
(2020). However, these works heavily rely on the distributional assumption of
the count data with the pre-specified dispersion parameter.

Most studies on NBR assumed the dispersion parameter as a constant. In
practice, however, not all models satisfy the assumption. Thus the need to
model the dispersion parameter as a function of some covariates. The
heterogeneous negative binomial regression (HNBR) extends the NBR by
observation-specific parameterization of the dispersion parameter (Hilbe,
2011). HNBR is a valuable tool for assessing the source of overdispersion. It
belongs to the double generalized linear models (DGLMs) or vector generalized
linear models (VGLMs), which are very useful in fitting more complex and
potentially realistic models (Yee, 2015; Nguelifack and Kemajou-Brown, 2019).
However, little work has been done to select the dispersion explanation
variables.

In this paper, we study the variable selection and dispersion estimation for
the heterogeneous NBR models. To the best of our knowledge and based on the
literature, this study is the first. Specifically, we propose a double
regression to estimate the coefficients of NB dispersion and NBR
simultaneously. Because of the high dimension of the covariates, we apply a
double $\ell_{1}$ penalty to both regressions. The two adjustment parameters
we set are different because the first-order conditions for estimating the
regression coefficients are entirely different from those for estimating the
dispersion parameters. We construct an algorithm to do variable selection and
dispersion estimation simultaneously. Similar studies on high-dimensional NBR
models include Wang et al. (2016), which assumed the dispersion parameter as a
constant. Their method requires an iterative algorithm to estimate the mean
regression and dispersion alternatively and implement lasso in each iteration.
If there are many iterations, such an algorithm is a waste of computing
resources.

The rest of the paper is organized as follows. Section 2 introduces the
heterogeneous overdispersed count data model and defines the double
$\ell_{1}$-penalized estimators for the mean and dispersion regressions. Then
we use a technique called the stochastic Lipschitz condition to derive the
asymptotic results in Section 3. Simulation studies and a real data
application are given in Section 4. Finally, section 5 concludes the article
with a discussion. All proofs and technical details are provided in the
Appendix.

##  2 Double $\ell_{1}$-penalized NBR

###  2.1 Heterogeneous overdispersed count data regressions

Suppose we have $n$ count responses ${Y_{i}}$ and $p$-dimensional covariates
$X_{i}=({x_{i1}},\cdots,{x_{ip}})$, $i\in[n]:=\\{1,2,\ldots,n\\}$. For the
Poisson regression models, the response obeys the Poisson distribution

$\mathrm{P}({Y_{i}=y_{i}}{\left|\lambda\right._{i}})=\frac{{\lambda_{i}^{{y_{i}}}}}{{{y_{i}}!}}{e^{-{\lambda_{i}}}},\qquad
i\in[n].$

With ${\lambda_{i}}={\rm{E}}{(Y_{i})}$, we require that the positive parameter
${\lambda_{i}}$ is related to a linear combination of $p$ covariates. A
plausible assumption for the link function is
$\eta({\lambda_{i}})=\log(\lambda_{i})={X_{i}^{\top}}\beta$. It is worth
noting that
$\mathrm{E}(Y_{i}|X_{i})=\mathrm{var}(Y_{i}|X_{i})=\exp({X_{i}^{\top}}\beta)>0.$

For the traditional negative binomial regression, it assumes that the count
data response obeys the NB distribution with over-dispersion:

$\mathrm{P}({Y_{\rm{i}}}={y_{i}}|{X}_{i})=:f({y_{i}};k,{\mu_{i}})=\frac{{\Gamma(k+{y_{i}})}}{{\Gamma(k){y_{i}}!}}{(\frac{{{\mu_{i}}}}{{k+{\mu_{i}}}})^{{y_{i}}}}{(\frac{k}{{k+{\mu_{i}}}})^{k}},\qquad
i\in[n],$ (2.1)

with $\mathrm{E}({Y_{i}}|X_{i})={\mu_{i}}=\exp(\beta^{\top}X_{i})$ and $k$ is
an unknown qualification of the overdispersion level. When $k\to\infty$, we
have
${\rm{var}}({Y_{i}}|X_{i})={\mu_{i}}+\frac{{\mu_{i}^{2}}}{k}\to{\mu_{i}}{\rm{=E(}}{Y_{i}}|X_{i})$,
the Poisson regression for the mean parameter ${\mu_{i}}$. Thus the Poisson
regression is a limiting case of negative binomial regression when the
dispersion parameter $k$ tends to infinite.

In the heterogeneous negative binomial regression, $k$ is proposed a specific
parameterization, i.e. $k=k(X_{i})$. More specifically, we assume in this
paper that

$\mu(x)=\exp\\{\theta^{(1)\top}x\\},\qquad k(x)=\exp\\{\theta^{(2)\top}x\\}.$

For notation simplicity, we denote

$Pf:=\mathrm{E}f(X_{i},Y_{i}),\qquad\mathbb{P}_{n}f:=\frac{1}{n}\sum_{i=1}^{n}f(X_{i},Y_{i}),\qquad\mathbb{G}_{n}f:=\sqrt{n}(\mathbb{P}_{n}-P)f,$

for any measurable function $f$.

Let $\theta=(\theta^{(1)\top},\theta^{(2)\top})^{\top}\in\mathbb{R}^{2p}$,the
log-likelihood is

$\displaystyle n\ell(\theta)$
$\displaystyle=\log\prod\limits_{i=1}^{n}{f({y_{i}},k_{i},{\mu_{i}})}=\sum\limits_{i=1}^{n}{\log}\left\\{\frac{{\Gamma(k_{i}+{Y_{i}})}}{{\Gamma(k_{i}){Y_{i}}!}}{(\frac{{{\mu_{i}}}}{{k_{i}+{\mu_{i}}}})^{{Y_{i}}}}{(\frac{k_{i}}{{k_{i}+{\mu_{i}}}})^{k_{i}}}\right\\}$
$\displaystyle=\sum\limits_{i=1}^{n}\bigg{[}-\log\frac{\Gamma\big{(}\exp\\{X_{i}^{\top}\theta^{(2)}\\}\big{)}}{\Gamma\big{(}{Y_{i}}+\exp\\{X_{i}^{\top}\theta^{(2)}\\}\big{)}}+{Y_{i}}X_{i}^{\top}\big{(}\theta^{(1)}-\theta^{(2)}\big{)}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad-\big{[}Y_{i}+\exp\\{X_{i}^{\top}\theta^{(2)}\\}\big{]}\log\Big{(}1+\exp\\{X_{i}^{\top}\big{(}\theta^{(1)}-\theta^{(2)}\big{)}\\}\Big{)}-\log{Y_{i}}!\bigg{]}$

We use the negative log-likelihood as the loss function $\gamma$, and define

$\displaystyle\gamma(\theta)$ $\displaystyle:=-\log f(y\,|\,x,\theta)+\log
y!.$

Denote $\partial_{j}:=\frac{\partial}{\partial\theta^{(j)}}$, $j=1,2$, the
score function for $\theta^{(1)}$ is

$\partial_{1}\ell(\theta)=\mathbb{P}_{n}\partial_{1}\gamma(\theta)=\frac{1}{n}\sum\limits_{i=1}^{n}{({Y_{i}}-{{\mathop{e}\nolimits}^{X_{i}^{\top}{\theta^{(1)}}}})\frac{{{{\mathop{e}\nolimits}^{X_{i}^{\top}{\theta^{(2)}}}}{X_{i}}}}{{{{\mathop{e}\nolimits}^{X_{i}^{\top}{\theta^{(1)}}}}+{{\mathop{e}\nolimits}^{X_{i}^{\top}{\theta^{(2)}}}}}}}.$

Furthermore, fix $\theta^{(1)}$, the score function for $\theta^{(2)}$ is

$\partial_{2}\ell(\theta)=\mathbb{P}_{n}\partial_{2}\gamma(\theta)=\frac{1}{n}\sum\limits_{i=1}^{n}{\left\\{{\left[{\log\left({1+{{\mathop{e}\nolimits}^{X_{i}^{\top}({\theta^{(1)}}-{\theta^{(2)}})}}}\right)-\sum\limits_{j=0}^{{Y_{i}}-1}{\frac{1}{{j+{{\mathop{e}\nolimits}^{X_{i}^{\top}{\theta^{(2)}}}}}}}}\right]{\rm{+}}\frac{{{Y_{i}}-{{\mathop{e}\nolimits}^{X_{i}^{\top}{\theta^{(1)}}}}}}{{{{\mathop{e}\nolimits}^{X_{i}^{\top}{\theta^{(1)}}}}+{{\mathop{e}\nolimits}^{X_{i}^{\top}{\theta^{(2)}}}}}}}\right\\}{{\mathop{e}\nolimits}^{X_{i}^{\top}{\theta^{(2)}}}}{X_{i}}}.$

It is easy to verify that

$P\partial_{1}\ell(\theta)=P\partial_{2}\ell(\theta)=0.$

So from now, we will suppose the true parameter is $\theta^{*}$.

###  2.2 Heterogeneous overdispersed NBR via double $\ell_{1}$-penalty

The lasso estimator under our circumstance is defined as

$\widehat{\theta}=\mathop{\operatorname{argmin}}_{\theta\in\Theta}\big{(}\mathbb{P}_{n}\gamma(\theta)+\lambda\|\theta\|_{\omega,1}\big{)},$
(2.2)

where $\lambda>0$ is the tuning parameter and the weighted norm is defined by

$\lambda\|\theta\|_{\omega,1}=\lambda_{1}\|\theta^{(1)}\|_{1}+\lambda_{2}\|\theta^{(2)}\|_{1}=\lambda\left(\omega_{1}\|\theta^{(1)}\|_{1}+\omega_{2}\|\theta^{(2)}\|_{1}\right),$

and
$\omega=(\omega_{1},\omega_{2})^{\top}=(\lambda_{1}/\lambda,\lambda_{2}/\lambda)^{\top}\in[0,1]\times[0,1]$
is the weight, $\|\cdot\|_{1}$ means the $\ell_{1}$-norm. This technique is
also used in Huang et al. (2021). Equation (2.2) is a weighted double
$\ell_{1}$-penalized problem, which is a kind of convex penalty optimization,
and when $\lambda_{1}=\lambda_{2}$, it becomes a single penalized problem. In
this paper, we use different $\lambda_{1}$ and $\lambda_{2}$, since the first-
order conditions for estimating the regression coefficients are entirely
different from those for estimating the dispersion parameters, and take
$\lambda=\lambda_{1}\vee\lambda_{2}$.

Since the weighted group lasso estimator $\widehat{\theta}_{n}$ has no closed-
form solution, we need to use iterative methods like quasi-Newton or
coordinate descent methods. We use BIC to choose the parameter $\lambda_{1}$
and $\lambda_{2}$.

$\operatorname{BIC}(\lambda_{1},\lambda_{2})=-2\ell(\widehat{\theta}_{n})+\dfrac{\log
n}{n}k,$

where $k$ is the number of nonzero estimated coefficients. To illustrate the
algorithm explicitly, we rewrite $\gamma(\theta)$ as
$\gamma(\theta^{(1)\top}x,\ \theta^{(2)\top}x)$, and define
$\theta^{(3)}=\lambda_{2}/\lambda_{1}\theta^{(2)}$,
$\theta^{{\dagger}}=(\theta^{(1)\top},\theta^{(3)\top})^{\top}$. Converting
$\theta^{(2)}$ into $\theta^{(3)}$ turns the double $\ell_{1}$-penalized
problem into a single penalized one, which is very to solve through some R
packages. The algorithm is formally given in Algorithm 1.

Algorithm 1 Double $\ell_{1}$-Penalized Optimization

Input: the set of tuning parameters
$\Lambda=\\{(\lambda_{1,i},\lambda_{2,i})\\}_{i=1}^{m}$
Output: the estimate $\widehat{\theta}_{n}$

for $i=1,...,m$, do

let $x^{*}=\frac{\lambda_{1,i}}{\lambda_{2,i}}x$;

solve
$\widehat{\theta}^{{\dagger}}=(\widehat{\theta}^{(1)\top},\widehat{\theta}^{(3)\top})^{\top}=\mathop{\operatorname{argmin}}_{\theta^{{\dagger}}\in\Theta}\big{(}\mathbb{P}_{n}\gamma(\theta^{(1)\top}x,\
\theta^{(3)\top}x^{*}))+\lambda_{1,i}\|\theta^{{\dagger}}\|_{1}\big{)}$;

obtain the estimate
$\widehat{\theta}_{n,i}=(\widehat{\theta}^{(1)\top},\frac{\lambda_{2,i}}{\lambda_{1,i}}\widehat{\theta}^{(3)\top})^{\top}$;

compute
$\operatorname{BIC}(\lambda_{1,i},\lambda_{2,i})=-2\ell(\widehat{\theta}_{n,i})+\dfrac{\log
n}{n}k_{i}$;

end for

find
$i_{opt}=\mathop{\operatorname{argmin}}_{i=1,...,m}{\operatorname{BIC}(\lambda_{1,i},\lambda_{2,i})}$;

return $\widehat{\theta}_{n,i_{opt}}$

The proposed algorithm can do variable selection and dispersion estimation
simultaneously. Similar studies on high-dimensional NBR models include Wang et
al. (2016), which assumed the dispersion parameter as a constant. However,
their method requires an iterative algorithm to estimate the mean regression
and dispersion alternatively and implement lasso in each iteration. If there
are many iterations, such an algorithm is a waste of computing resources.

##  3 Main results

###  3.1 Stochastic Lipschitz conditions

We write the maximum of $Y_{i}$ from the sample of size $n$ as $M_{Y,n}$, then
the sample space for $\\{Y_{i}\\}_{i=1}^{n}$ is
$\mathcal{Y}:=\\{y\in\mathbb{N},\,y\leq M_{y,n}\\}$. i.e.
$M_{y,n}=\max_{i\in[n]}Y_{i}$. Note that
$\lim_{n\rightarrow\infty}\mathrm{P}(M_{y,n}=\infty)=1$, what we need to
tackle is actually a unbounded empirical process. But for
$z:=\left(\begin{matrix}x^{\top}&\\\
&x^{\top}\end{matrix}\right)\in\mathbb{R}^{2\times 2p}$, we can assume the
value space $\mathcal{S}$ for $s:=z\theta$ is bounded and satisfies

$\mathcal{S}:=\big{\\{}s=(s_{1},s_{2})^{\top}\in\mathbb{R}^{2},\
-\infty<m_{s,n}\leq s_{j}\leq|s_{j}|\leq M_{s,n}<\infty,\ j=1,2\big{\\}}.$

As we can see, the most significant difference between this article and other
conventional literature about lasso estimators is that we use $s=z\theta$
rather than $\theta$ as the explanatory variable to analyze the properties of
the loss function $\gamma$. This is not a traditional way. At the first
glimpse, the combination may complicate the analysis in the next step because
the KKT condition requires the story about
$\frac{\partial}{\partial\theta}\gamma$, which is critical for the traditional
convex penalty problem. However, this article will try a different approach,
the stochastic Lipschitz conditions introduced in the event $\mathcal{A}$ of
Proposition 1 in Zhang and Jia (2022), to solve $\ell_{1}$-penalization
problem. Define the _local stochastic Lipschitz constant_ by

$\operatorname{Lip}(f;\theta^{*}):=\sup_{\theta\in\Theta/\\{\theta^{*}\\}}\left|\frac{\sqrt{n}\mathbb{G}_{n}\big{(}f(\theta)-f(\theta^{*})\big{)}}{\|\theta-\theta^{*}\|_{1}}\right|.$

The most apparent advantage of the stochastic Lipschitz conditions over KKT
condition is that it can easily deal with the several parameters involved in
different locations of the model that need to impose the same penalty on them,
which is why we do not need to derive KKT condition in this paper.

To establish the stochastic Lipschitz conditions for this unbounded counting
process, another assumption called strongly midpoint log-convex for some
positive $\gamma$ should be satisfied, which states for the joint density from
the sample $\mathbb{Y}:=(Y_{1},\ldots,Y_{n})^{\top}\in\mathbb{Z}^{n}$’s
negative log-density of $n$ independent NB responses $\psi({y}):=-\log
p_{\mathbb{Y}}({y})$ satisfies

$\psi({x})+\psi({y})-\psi\left(\left\lceil\frac{1}{2}{x}+\frac{1}{2}{y}\right\rceil\right)-\psi\left(\left\lfloor\frac{1}{2}{x}+\frac{1}{2}{y}\right\rfloor\right)\geq\frac{\gamma}{4}\|{x}-{y}\|_{2}^{2},\qquad\forall{x},{y}\in\mathbb{Z}^{n}.$

This assumption is a condition that ensures that the suprema of the multiplier
empirical processes of $n$ independent responses have sub-exponential
concentration phenomena, which can be alternatively be checked by the tail
inequality for suprema of empirical processes corresponding to classes of
unbounded functions (Adamczak (2008)).

######  Corollary 3.1.

Suppose $\max_{i\in[n],\,1\leq k\leq p}|X_{ik}|\leq M_{x}<\infty$, the
parameter space $\Theta$ is convex and its diameter $D_{\Theta}<\infty$. If
$\\{Y_{i}\\}_{i=1}^{n}$ and $\\{Z_{i}\theta\\}_{i\in[n],\,\theta\in\Theta}$
are both in the value space $\mathcal{Y}$ and $\mathcal{S}$ defined as
previous, then for any $\theta\in\Theta$,

$\displaystyle\operatorname{Lip}(\gamma;\theta^{*})$
$\displaystyle=\sup_{\theta\in\Theta/\\{\theta^{*}\\}}\left|\frac{\sqrt{n}\mathbb{G}_{n}\big{(}\gamma(\theta)-\gamma(\theta^{*})\big{)}}{\|\theta-\theta^{*}\|_{1}}\right|$
$\displaystyle\leq\sqrt{n}M_{q}:=\Big{(}A_{1}\sqrt{\log(2p/q_{2})}+A_{2}\sqrt{\log
p}+A_{3}\sqrt{\log(p/q_{3})}\Big{)}\sqrt{\max_{1\leq k\leq
p}\sum_{i=1}^{n}X_{ik}^{2}}$
$\displaystyle+B\sqrt{\log(2p/q_{1})}\sqrt{\Big{(}\max_{1\leq k\leq
p}\sum_{i=1}^{n}X_{ik}^{4}\Big{)}^{1/2}}\vee C\log(2p/q_{1})+D\log(p/q_{3}),$

with probability at least $1-q_{0}$, where $q_{1},q_{2},q_{3}\in(0,1)$ satisfy
$q_{1}+q_{2}+q_{3}=q_{0}$, and the constants are as follows:

$\displaystyle A_{1}=\sqrt{2}F_{1},\qquad
A_{2}=32\sqrt{2}M_{x}F_{2}D_{\Theta},\qquad
A_{3}=\sqrt{2}\big{(}2(F_{1}+M_{y,n})\vee F_{2}M_{x}D_{\Theta}\big{)},$
$\displaystyle B=6\sqrt{2\big{(}w^{(1)}\vee
w^{(2)}\big{)}\bigg{(}\sum_{i=1}^{n}a(\mu_{i},k_{i})^{4}\bigg{)}^{1/2}},\qquad
C=12M_{x}\big{(}w^{(1)}\vee w^{(2)}\big{)}\max_{1\leq i\leq
n}a(\mu_{i},k_{i}),$ $\displaystyle D=8\big{(}2(F_{1}+M_{y,n})\vee
F_{2}M_{x}D_{\Theta}\big{)}M_{x},\qquad
w^{(1)}=\frac{e^{M_{s,n}}}{e^{m_{s,n}}+e^{M_{s,n}}},~{}~{}w^{(2)}=\frac{e+e^{M_{s,n}-m_{s,n}}}{1+e^{m_{s,n}-M_{s,n}}}+\frac{1}{1+e^{m_{s,n}-M_{s,n}}},$

where $M_{y,n}=\max_{i\in[n]}Y_{i}$ is the suprema empirical process.

It is worthy to note that the $M_{y,n}$ in Corollary 3.1 is a random process,
hence the bound above is not deterministic. Fortunately, $M_{y,n}$ can use the
strongly midpoint log-convex condition to be bounded, which we state in Lemma
A.3. Corollary 3.1 combined Lemma A.3 will give the following corollary as a
step more.

######  Corollary 3.2.

Assume the conditions are the same as that in Corollary 3.1, then the
stochastic Lipschitz constant has a nonrandom upper bound:

$\displaystyle\operatorname{Lip}(\gamma;\theta^{*})$
$\displaystyle\leq\sqrt{n}M_{q}^{\prime}:=\Big{(}A_{1}\sqrt{\log(2p/q_{2})}+A_{2}\sqrt{\log
p}+2A_{3}^{\prime}\big{(}\log(2n/q_{4})+\sqrt{\log(np/q_{3})\big{)}}\Big{)}\sqrt{\max_{1\leq
k\leq p}\sum_{i=1}^{n}X_{ik}^{2}}$
$\displaystyle+B\sqrt{\log(2p/q_{1})}\sqrt{\Big{(}\max_{1\leq k\leq
p}\sum_{i=1}^{n}X_{ik}^{4}\Big{)}^{1/2}}\vee C\log(2p/q_{1})+D\log(p/q_{3}),$

with probability at least $1-q_{0}$, where $q_{1},q_{2},q_{3},q_{4}\in(0,1)$
satisfy $q_{1}+q_{2}+q_{3}+q_{4}=q_{0}$, and

$A_{3}^{\prime}=2\sqrt{2}\left(F_{1}+\left(2\gamma\max_{i\in[n]}\left[a(\mu_{i},k_{i})-\frac{\mu_{i}}{\log
2}\right]\right)\vee F_{2}M_{x}D_{\Theta}\right).$

Theorem 3.1 gives us a different sight of the loss function far more than KKT
conditions. However, the stochastic Lipschitz condition above does not compare
the estimated and true values directly. We can resolve this issue by using an
eigenvalue condition on the design matrix consisting of $X_{i}$. Since the
design matrix $\mathrm{X}$ is fixed, the eigenvalue condition in the next
section is reasonable. It is worthy to note that this inequality is an oracle
since it involves an unknown empirical process on the right side.

###  3.2 $\ell_{2}$-estimation error oracle inequalities RE conditions

As we said previously, although we use stochastic Lipschitz conditions instead
of KKT conditions, the restricted eigenvalue conditions (RE conditions) are
still required. We denote by $\delta_{J}$ the vector in $\mathbb{R}^{p}$ with
the same coordinates as $v$ on $J$ and zero coordinates on the complement
$J^{c}$ of $J$, and $\operatorname{spt}(v)=\\{j:v_{j}\neq 0\\}$. We will
assume that the minima in (2.2) can always be obtained in the following
setting, but it may not be unique. In general, to bound
$\widehat{\theta}-\theta^{*}$, some conditions on the design matrix
$\mathrm{X}\in\mathbb{R}^{n\times p}$ are needed for getting abound in terms
of the $\ell_{2}$ norm of $\theta-\theta^{*}$. Here we will utilize the
restricted eigenvalue condition introduced in Bickel et al. (2009), which says
that for some $1\leq s\leq p$ and $K>0$,

$\kappa(s,K)=\min\left\\{\frac{\|\mathrm{X}v\|_{2}}{\sqrt{n}\|v_{J}\|_{2}}:1\leq|J|\leq
s,\,v\in\mathbb{R}^{p}/\\{0\\},\,\|v_{J^{c}}\|_{1}\leq
K\|v_{J}\|_{1}\right\\}>0.$ (3.1)

It should be noted that omitting the weight $\omega$ and the sparse restricted
set $\|v_{J^{c}}\|_{1}\leq K\|v_{J}\|_{1}$ leads to
$v^{\top}\big{[}\frac{1}{n}\mathrm{X}^{\top}\mathrm{X}\big{]}v/v^{\top}v\geq\kappa^{2}(s,K)$.
Thus it means that the smallest eigenvalue of the sample covariance matrix
$\frac{1}{n}\mathrm{X}^{\top}\mathrm{X}$ is positive, which is impossible when
$p>n$ since $\frac{1}{n}\mathrm{X}^{\top}\mathrm{X}$ is not full rank. To
avoid this problem, Bickel et al. (2009) consider the restricted eigenvalue
condition under the sparse restricted set $\|v_{J^{c}}\|_{1}\leq
K\|v_{J}\|_{1}$ as considerable relation in sparse high-dimensional
estimation. The restricted eigenvalue is from the restricted strong convexity,
which enforces a strong convexity condition for the negative log-likelihood
function of linear models under certain sparse restrict set.

Due to the double penalty, besides the RE condition, we also require another
condition similar to RE condition so-called $l$-restricted isometry constant
defined in Candes et al. (2007) as follows

$\sigma_{\mathrm{X},l}^{2}=\max\left\\{{\|\mathrm{X}v\|_{2}^{2}}/{\|v\|_{2}^{2}}:v\in\mathbb{R}^{p},\,1\leq\operatorname{spt}(v)\leq
l\right\\}\in(0,\infty),$

which essentially requires the eigenvalue of sample covariance matrix under
every vector with cardinality less than $l$ ($l$ should be no more than $n$)
approximately behaves normally like the low dimensional case.

With the RE condition and $l$-restricted isometry constant, and the two
theorems we established before, the lasso estimator in (2.2) can guarantee a
good consistent property.

######  Lemma 3.1 (see Lemma 3.1 in Candes et al. (2007)).

Suppose $T_{0}$ is a set of cardinality $S$. For a vector
$h\in\mathbb{R}^{p}$, we let $T_{1}$ be the $S$ largest positions of $h$
outside of $T_{0}$. Put $T_{01}=T_{0}\cup T_{1}$, then

$\|h\|_{2}^{2}\leq\|h_{T_{01}}\|_{2}^{2}+S^{-1}\|h_{T_{0}^{c}}\|_{1}^{2}.$

######  Theorem 3.1.

Suppose the condition is the same as that in Theorem A.1. Furthermore, assume
$p_{1}=\operatorname{spt}\big{(}\theta^{(*1)}\big{)}\vee\operatorname{spt}\big{(}\theta^{(*2)}\big{)}\leq
p/2$, and there exists some $K>1$, $\kappa:=\kappa(2p_{1},K)>0$. Let
$\lambda=\frac{(K+1)M_{q}}{n(K-1)}$, then using this $\lambda$ in (2.2), with
probability at least $1-q$,

$\|\widehat{\theta}-\theta^{*}\|_{2}^{2}\leq\frac{8p_{1}M_{q}^{2\prime}K^{2}}{\kappa^{4}n^{2}C_{\gamma}^{2}(K-1)^{2}}\left[2+K^{2}+\frac{2(1+2p_{1}K^{2})(n\kappa^{2}+2\sigma_{\mathrm{X},p_{1}}^{2})}{n\kappa^{2}}\right],$

where $M_{q}$, $C_{\gamma}$ are defined in Theorem 3.1 and A.1 respectively.

######  Remark 3.1.

Compared to the single lasso problem, in which we only have one unknown
vectorized parameter, the oracle inequality in Theorem 3.1 has an extra term
$\frac{2(1+2p_{1}K^{2})(n\kappa^{2}+2\sigma_{\mathrm{X},p_{1}}^{2})}{n\kappa^{2}}$.

######  Remark 3.2.

From Theorem 3.1, we know that the $\ell_{2}$ convergence rate is minimax
optimal, as studied in Zhang and Jia (2022).

From Theorem 3.1, we can get the following corollary regarding the consistency
of $\widehat{\theta}$ immediately by the property of the supreme empirical
process.

##  4 Numerical Studies

###  4.1 Simulations

In this section, we evaluate the finite sample performance of the proposed
method. The response is generated from the negative binomial regression model
(2.1) with

$\mu(x)=\exp\\{\theta^{(1)\top}x\\},\text{and~{}}k(x)=\exp\\{\theta^{(2)\top}x\\},$

where $\theta^{(1)}$ and $\theta^{(2)}$ are two $p$-dimensional parameters.
The explanatory variables are generated from the multivariate normal
distributions with mean vector $\bf{0}$ and $Cov(x_{i},x_{j})=\rho^{|i-j|}$,
where $\rho=0,0.5$. The following two examples show the performance of the
proposed estimator for the low dimensional heterogeneous negative binomial
regression and the variable selection in the high dimensional case,
respectively. The R package lbfgs is required to solve the optimization
problem.

Example 1 (Low dimension). We set $p=3$ and $n=100,200,400$. The true
parameters are $\theta^{(1)}=(1,2,-1)$ and $\theta^{(2)}=(-1,0.5,1)$, and
their maximum likelihood estimators are denoted as $\hat{\theta}^{(1)}$ and
$\hat{\theta}^{(2)}$, respectively. We compare the estimator
$\hat{\theta}^{(1)}$ with $\hat{\theta}^{(1)*}$, which ignores the
heterogeneity of the overdispersion and treats $k(x)$ as a constant. Table 1
displays the average squared estimation errors
$\|\hat{\theta}-{\theta}\|_{2}^{2}$ based on 200 repetitions.

We can make the following observations from the table. Firstly, the
performances of the three estimators become better and better as $n$
increases. Secondly, the estimator $\hat{\theta}^{(1)}$, which estimates the
parameter in the mean function $\mu(x)$, performs better than
$\hat{\theta}^{(2)}$, which estimates the parameter in the overdispersion
function $k(x)$. Last but the most important, $\hat{\theta}^{(1)*}$ performs
much worse than $\hat{\theta}^{(1)}$. For example, the average squared
estimation error of $\hat{\theta}^{(1)*}$ is about 5 times of
$\hat{\theta}^{(1)}$’s when $n=100$, and 10 times of $\hat{\theta}^{(1)}$’s
when $n=400$. The comparison between $\hat{\theta}^{(1)}$ and
$\hat{\theta}^{(1)*}$ indicates the necessity of considering the heterogeneity
of the overdispersion.

Table 1:  The average squared estimation errors of the estimators. n | $\rho=0$ | $\rho=0.5$
---|---|---
| $\hat{\theta}^{(1)*}$ | $\hat{\theta}^{(1)}$ | $\hat{\theta}^{(2)}$ | $\hat{\theta}^{(1)*}$ | $\hat{\theta}^{(1)}$ | $\hat{\theta}^{(2)}$
100 | 0.1597 | 0.0335 | 0.72414 | 0.1809 | 0.0397 | 0.68904
200 | 0.0862 | 0.01 | 0.22149 | 0.0837 | 0.0169 | 0.33048
400 | 0.05 | 0.0047 | 0.08847 | 0.0619 | 0.0067 | 0.15066

Example 2 (High dimension). The sample sizes are chosen to be $n=100,200,400$,
with dimension $p\in(25,50,150)$, $(50,100,250)$ and $(100,200,500)$,
respectively. We set $\theta^{(1)}=(1,2,-1,0,\dots,0)$ and
$\theta^{(2)}=(-1,0.5,1,0,\dots,0)$. The unknown tuning parameters
$(\lambda_{1},\lambda_{2})$ for the penalty functions are chosen by BIC
criterion in the simulation. Results over 200 repetitions are reported. For
each case, table 2 reports the number of repetitions that each important
explanatory variable is selected in the final model and also the average
number of unimportant explanatory variables being selected.

It shows that the variable selection procedure performs better and better as
the sample size $n$ increases. When $n=400$, the important explanatory
variables in $\mu(x)$ and $k(x)$ are correctly selected in almost every
repetitions. When the dimension $p$ increases, the procedure may select more
unimportant explanatory variables, but the average numbers are less than
$1.3$. The important variables in $k(x)$ are less likely to be selected than
the important variables in $\mu(x)$ especially when the sample size is small,
as well as the unimportant variables.

Table 2:  The results of variable selection. |  | $\mu(x)$ | $k(x)$
---|---|---|---
n | p | $\theta^{(1)}_{1}$ | $\theta^{(1)}_{2}$ | $\theta^{(1)}_{3}$ | Other $\theta^{(1)}$s | $\theta^{(2)}_{1}$ | $\theta^{(2)}_{2}$ | $\theta^{(2)}_{3}$ | Other $\theta^{(2)}$s
$\rho=0$
100 | 25 | 192 | 200 | 190 | 0.37 | 180 | 184 | 180 | 0.32
| 50 | 196 | 200 | 193 | 0.52 | 182 | 180 | 188 | 0.41
| 150 | 194 | 194 | 192 | 1.02 | 188 | 182 | 186 | 0.41
200 | 50 | 200 | 200 | 200 | 0.59 | 200 | 190 | 198 | 0.53
| 100 | 200 | 200 | 200 | 0.91 | 196 | 186 | 196 | 0.69
| 250 | 199 | 199 | 198 | 1.18 | 198 | 198 | 198 | 0.69
400 | 100 | 200 | 200 | 200 | 0.4 | 200 | 198 | 200 | 0.55
| 200 | 200 | 200 | 200 | 0.6 | 200 | 200 | 200 | 0.51
| 500 | 200 | 200 | 200 | 1.21 | 200 | 200 | 200 | 0.61
$\rho=0.5$
100 | 25 | 194 | 198 | 194 | 0.41 | 179 | 184 | 180 | 0.35
| 50 | 196 | 196 | 190 | 0.63 | 178 | 182 | 180 | 0.42
| 150 | 194 | 196 | 192 | 1.01 | 180 | 184 | 182 | 0.43
200 | 50 | 200 | 200 | 198 | 0.38 | 196 | 183 | 190 | 0.32
| 100 | 199 | 200 | 198 | 0.53 | 194 | 186 | 194 | 0.44
| 250 | 196 | 198 | 196 | 1.1 | 196 | 196 | 194 | 0.55
400 | 100 | 200 | 200 | 200 | 0.28 | 200 | 199 | 194 | 0.34
| 200 | 200 | 200 | 198 | 0.47 | 200 | 198 | 196 | 0.36
| 500 | 200 | 200 | 198 | 1.07 | 200 | 198 | 196 | 0.56

###  4.2 A real data example

In this section, we apply the proposed method to the dataset of German health
care demand. The data was employed in Riphahn et al. (2003) and could be
downloaded on http://qed.econ.queensu.ca/jae/2003-v18.4/riphahn-wambach-
million/. The data contains 27,326 observations on 25 variables, including two
dependent variables Docvis (number of doctor visits in last three months) and
Hospvis (number of hospital visits in last calendar year). For conciseness, we
focus on Docvis in this study. We build the HNBR model based on the proposed
variable selection procedure and make the standard NBR model a comparison.
Define the fitting errors (FE) as
$n^{-1}\sum_{i=1}^{n}(y_{i}-\widehat{y}_{i})$, where $y_{i}$ denotes the raw
data of Hospvis, $\hat{y}_{i}$ is the predicted value, and $n$ is the sample
size. As the data is observed during 1984-1988, 1991, and 1994, we make the
analysis each observed year. Table 3 displays the variable selection results
and fitting errors.

We have the following findings from the table. First, the important variables
in the NBR are the same as HNBR models in each year, and the estimates are
close. Second, the selected variables in $\mu(x)$ are almost the same every
year, namely Age, Hsat (health satisfaction), Handper (degree of handicap),
and Educ (years of schooling). Moreover, some of these variables still play an
essential role in $k(x)$, and $k(x)$ contains no variables other than these.
Also we can see that the fitting errors of HNBR is less than that of NBR. All
of these illustrate the advantage of our method.

Table 3:  The variable selection results and the fitting errors (FE) of NBR and HNBR models. The variable Others $=\\{$ Married, Haupts, Reals, Fachhs, Abitur, Univ, Working, Bluec, Whitec, Self, Beamt, Public, Addon $\\}$. Since these variables are not selected in any year, we put them in “Others” for brevity. Variables | $1984$ | 1985 | $1986$ | 1987
---|---|---|---|---
| NBR | HNBR | NBR | HNBR | NBR | HNBR | NBR | HNBR
|  | $\mu(x)$ | $k(x)$ |  | $\mu(x)$ | $k(x)$ |  | $\mu(x)$ | $k(x)$ |  | $\mu(x)$ | $k(x)$
Female | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0
Age | -0.013 | -0.013 | -0.012 | -0.009 | -0.01 | -0.007 | -0.006 | -0.006 | -0.013 | -0.002 | -0.001 | -0.018
Hsat | -0.205 | -0.2 | -0.025 | -0.244 | -0.237 | 0 | -0.188 | -0.195 | -0.045 | -0.158 | -0.153 | -0.043
Handdum | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0
Handper | 0.005 | 0.005 | 0.004 | 0.007 | 0.006 | 0.007 | 0.007 | 0.007 | 0 | 0.007 | 0.007 | 0.01
Hhninc | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0
Hhkids | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0
Educ | 0 | 0 | -0.027 | 0 | 0 | -0.064 | -0.035 | -0.038 | 0 | -0.095 | -0.106 | -0.003
Others | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0
FE | 0.798 | 0.602 | 2.203 | 1.874 | 0.735 | 0.581 | 1.314 | 1.027
Variables | $1988$ | 1991 | $1994$ |
| NBR | HNBR | NBR | HNBR | NBR | HNBR |  |  |
|  | $\mu(x)$ | $k(x)$ |  | $\mu(x)$ | $k(x)$ |  | $\mu(x)$ | $k(x)$ |  |  |
Female | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |  |  |
Age | -0.015 | -0.014 | -0.012 | -0.022 | -0.019 | -0.003 | -0.005 | -0.004 | -0.011 |  |  |
Hsat | -0.191 | -0.187 | -0.015 | -0.112 | -0.132 | -0.049 | -0.226 | -0.224 | -0.06 |  |  |
Handdum | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |  |  |
Handper | 0.011 | 0.009 | 0.006 | 0.014 | 0.013 | 0 | 0.007 | 0.008 | 0.004 |  |  |
Hhninc | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |  |  |
Hhkids | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |  |  |
Educ | -0.016 | -0.023 | -0.002 | -0.074 | -0.068 | 0 | -0.064 | -0.069 | 0 |  |  |
Others | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |  |  |
FE | 1.144 | 0.912 | 1.007 | 0.787 | 0.713 | 0.58 |  |  |

##  5 Conclusion and future study

We study the high-dimensional heterogeneous overdispersed count data via
negative binomial regression models, and propose a double
$\ell_{1}$-regularized method for simultaneous variable selection and
dispersion estimation. Under the restricted eigenvalue conditions, we prove
the oracle inequalities with lasso estimators of two partial regression
coefficients for the first time, using concentration inequalities of empirical
processes. Furthermore, we derive the consistency and convergence rate for the
estimators, which are the theoretical guarantees for further statistical
inference. Simulation studies and a real example from the German health care
demand data indicate that the proposed method works satisfactorily.

In this work, we assume that the responses are independent. In the time-series
data, however, the NB responses are temporal dependent (Yang et al., 2021).
Thus, weak dependence conditions, including $\rho$-mixing, $m$-dependent
types, could be considered in the future. The consequences of the main result
would motivate further study for statistical inference, such as testing
heterogeneous

$H_{0}:\theta^{(2)}=0\text{ vs. }H_{1}:\theta^{(2)}\neq 0.$

The issues concerning the hypothesis testing are via the debiased Lasso
estimator; see Shi et al. (2019) and references therein. Another possible
study is the false discovery rate (FDR) control, which aims to identify some
small number of statistically significantly nonzero results after getting the
sparse penalized estimation of HNBR; see Xie and Lederer (2021); Cui et al.
(2021).

##  Appendix A Proofs

The first step is giving the property of the loss function. From mathematical
analysis, we prefer bounded things to unlimited things. Denote $\partial_{j}$
is the first partial differentiation with respect to $s_{j}$. The bounded
aspect for $y$ and $s$ gives a nice property for the loss function
$\gamma(s,y)=\gamma(z\theta,y)$.

######  Lemma A.1.

We have

$\partial_{1}\gamma(s,y)=-\frac{ye^{s_{2}}}{e^{s_{1}}+e^{s_{2}}}+\frac{e^{s_{1}+s_{2}}}{e^{s_{1}}+e^{s_{2}}},\qquad\partial_{2}\gamma(s,y)=\nu(s,y)+\frac{ye^{s_{2}}}{e^{s_{1}}+e^{s_{2}}}$

where
$\nu(s,y)=e^{s_{2}}\big{(}\psi(e^{s_{2}})-\psi(y+e^{s_{2}})\big{)}+e^{s_{2}}\log\big{(}1+e^{s_{1}-s_{2}}\big{)}-\frac{e^{s_{1}+s_{2}}}{e^{s_{1}}+e^{s_{2}}}$
satisfying

$\sup_{s\in\mathcal{S},\,y\in\mathcal{Y}}|\nu(s,y)|\leq
F_{1},\qquad\sup_{s\neq
t\in\mathcal{S},\,y\in\mathcal{Y}}\frac{|\nu(s,y)-\nu(t,y)|}{\|s-t\|_{\infty}}\leq
F_{2},$

with
$F_{1}=M_{y,n}\big{(}1+e^{-m_{s,n}}\big{)}+e^{M_{s,n}}+\frac{e^{2M_{s,n}}}{2e^{m_{s,n}}}$,
and

$\displaystyle F_{2}$
$\displaystyle=2\left[\left|e^{M_{s,n}}\left(1+\log(M_{y,n}+e^{M_{s,n}})-\frac{1}{2(M_{y,n}+e^{M_{s,n}})}\right)\right|\vee\left|1-m_{s,n}e^{m_{s,n}}\right|\right]$
$\displaystyle\qquad\qquad\qquad+\left[\frac{e^{2M_{s,n}}}{e^{m_{s,n}}}+\frac{2e^{2M_{s,n}}}{e^{m_{s,n}}+e^{M_{s,n}}}\right]+\frac{3}{2}e^{M_{s,n}}.$

###### Proof.

We will use the properties of the psi function, the logarithmic derivative of
the gamma function, to prove this lemma. Write
$\psi(x)=\Gamma^{\prime}(x)/\Gamma(x)$. For any $s\in\mathcal{S},\
y\in\mathcal{Y}$, using the Binet ’s formula (see p.18 of Bateman (1953))

$\psi(x)=\log x-\int_{0}^{\infty}\varphi(t)e^{-tx}\,dt,$

where $\varphi(t)=1/(1-e^{-t})-1/t$ is strictly increasing on $(0,\infty)$, it
gives

$0<\psi^{\prime}(x)=\frac{1}{x}+\int_{0}^{\infty}t\varphi(t)e^{-tx}\,dt\leq\frac{1}{x}+\int_{0}^{\infty}te^{-tx}\,dt=\frac{1}{x}+\frac{1}{x^{2}}.$

and $y\geq 0$, we have

$\displaystyle|\nu(s,y)|$
$\displaystyle=\left|e^{s_{2}}\big{(}\psi(e^{s_{2}})-\psi(y+e^{s_{2}})\big{)}+e^{s_{2}}\log\big{(}1+e^{s_{1}-s_{2}}\big{)}-\frac{e^{s_{1}+s_{2}}}{e^{s_{1}}+e^{s_{2}}}\right|$
$\displaystyle\leq
e^{s_{2}}y\left(\frac{1}{e^{s_{2}}}+\frac{1}{e^{s_{2}}}\right)+e^{s_{2}}e^{s_{1}-s_{2}}+\frac{e^{s_{1}+s_{2}}}{e^{s_{1}}+e^{s_{2}}}\leq
M_{y,n}\big{(}1+e^{-m_{s,n}}\big{)}+e^{M_{s,n}}+\frac{e^{2M_{s,n}}}{2e^{m_{s,n}}}.$

Then the first inequality in the lemma has been verified. On the other hand,
by using the fact that (see (2.2) in Alzer (1997))

$\frac{1}{2x}<\log(x)-\psi(x)<\frac{1}{x},\qquad x>0,$

for the function $f_{1}(x)=e^{x}\psi(e^{x})$ and
$f_{2}(x)=e^{x}\psi(y+e^{x})$,

$\displaystyle
f_{1}^{\prime}(x)=e^{x}\psi(e^{x})+e^{2x}\psi^{\prime}(e^{x})\leq
e^{x}\left(x-\frac{1}{2e^{x}}\right)+e^{2x}\left(\frac{1}{e^{x}}+\frac{1}{e^{2x}}\right)=(x+1)e^{x}+\frac{1}{2},$
$\displaystyle
f_{1}^{\prime}(x)=e^{x}\psi(e^{x})-e^{2x}\psi^{\prime}(e^{x})\geq
e^{x}\psi(e^{x})\geq xe^{x}-1,$ $\displaystyle
f_{2}^{\prime}(x)=e^{x}\psi(y+e^{x})+e^{2x}\psi^{\prime}(y+e^{x})\leq
e^{x}\left(1+\log(y+e^{x})-\frac{1}{2(y+e^{x})}\right)+1,$ $\displaystyle
f_{2}^{\prime}(x)\geq e^{x}\left(\log(y+e^{x})-\frac{1}{y+e^{x}}\right)\geq
xe^{x}-1,$

and for any $s\neq t\in\mathcal{S}$, $y\in\mathcal{Y}$, we conclude that

$|\psi(e^{s_{2}})e^{s_{2}}-\psi(e^{t_{2}})e^{t_{2}}|=|f_{1}(s_{2})-f_{1}(t_{2})|\leq\left(\left|(M_{s,n}+1)e^{M_{s,n}}+{1}/{2}\right|\vee\left|1-m_{s,n}e^{m_{s,n}}\right|\right)\|s-t\|_{\infty},$

and

$\displaystyle|\psi(y+e^{s_{2}})e^{s_{2}}-\psi(y+e^{t_{2}})e^{t_{2}}|=|f_{2}(s_{2})-f_{2}(t_{2})|$
$\displaystyle\leq\left[\left|e^{M_{s,n}}\left(1+\log(M_{y,n}+e^{M_{s,n}})-\frac{1}{2(M_{y,n}+e^{M_{s,n}})}\right)\right|\vee\left|1-m_{s,n}e^{m_{s,n}}\right|\right]\|s-t\|_{\infty}.$

Besides, using the median value theorem again, we also have

$\displaystyle\left|e^{s_{2}}\log(1+e^{s_{1}-s_{2}})-e^{t_{2}}\log(1+e^{t_{1}-t_{2}})\right|$
$\displaystyle\leq$
$\displaystyle\log(1+e^{s_{1}-s_{2}})|e^{s_{2}}-e^{t_{2}}|+e^{t_{2}}\left|\log(1+e^{s_{1}-s_{2}})-\log(1+e^{t_{1}-t_{2}})\right|$
$\displaystyle\leq$ $\displaystyle
e^{2M_{s,n}-m_{s,n}}|s_{2}-t_{2}|+e^{M_{s,n}}\frac{1}{1+e^{-(M_{s,n}-m_{s,n})}}|(s_{1}-s_{2})-(t_{1}-t_{2})|$
$\displaystyle\leq$
$\displaystyle\left[\frac{e^{2M_{s,n}}}{e^{m_{s,n}}}+\frac{2e^{2M_{s,n}}}{e^{m_{s,n}}+e^{M_{s,n}}}\right]\|s-t\|_{\infty},$

and

$\displaystyle\left|\frac{e^{s_{1}+s_{2}}}{e^{s_{1}}+e^{s_{2}}}-\frac{e^{t_{1}+t_{2}}}{e^{t_{1}}+e^{t_{2}}}\right|$
$\displaystyle\leq
e^{s_{1}}\left|\frac{1}{1+e^{s_{1}-s_{2}}}-\frac{1}{1+e^{t_{1}-t_{2}}}\right|+\frac{1}{1+e^{t_{1}-t_{2}}}|e^{s_{1}}-e^{t_{1}}|$
$\displaystyle\leq e^{M_{s,n}}\frac{1}{4}|(s_{1}-s_{2})-(t_{1}-t_{2})|+1\times
e^{M_{s,n}}|s_{1}-t_{1}|\leq\frac{3}{2}e^{M_{s,n}}\|s-t\|_{\infty},$

where the fact used is that $f_{3}(x)=1/(1+e^{x})$ satisfies
$|f_{3}^{\prime}(x)|=1/(e^{x}+e^{-x}+2)\leq 1/4$. Since

$\displaystyle|\partial_{2}\gamma(s,y)-\partial_{2}\gamma(t,y)|$
$\displaystyle\leq\left|\psi(e^{s_{2}})e^{s_{2}}-\psi(e^{t_{2}})e^{t_{2}}\right|+\left|\psi(y+e^{s_{2}})e^{s_{2}}-\psi(y+e^{t_{2}})e^{t_{2}}\right|$
$\displaystyle\qquad+\left|e^{s_{2}}\log(1+e^{s_{1}-s_{2}})-e^{t_{2}}\log(1+e^{t_{1}-t_{2}})\right|+\left|\frac{e^{s_{1}+s_{2}}}{e^{s_{1}}+e^{s_{2}}}-\frac{e^{t_{1}+t_{2}}}{e^{t_{1}}+e^{t_{2}}}\right|,$

we can conclude the second inequality in the lemma. ∎

The Lemma separates the partial derivative of $\gamma$ into two parts: the
first part is the linear about the response variable $y$ (say
$-ye^{s_{2}}/(e^{s_{1}}+e^{s_{2}})$, $e^{s_{1}+s_{2}}/(e^{s_{1}}+e^{s_{2}})$,
and $ye^{s_{2}}/(e^{s_{1}}+e^{s_{2}})$), the second part is other complicated
functions (not linear function) about $y$. The first part is relatively easy
to analyze since the following concentration inequality gives a measure of
dispersion about the weighted summation of negative binomial variables. This
concentration inequality is a special case for the weighted summation of a
series of random variables, which can be proved by sub-exponential
concentration results in Proposition 4.2 in Zhang and Chen (2021).

######  Lemma A.2.

Suppose $\\{Y_{i}\\}_{i=1}^{n}$ are independently distributed as
$\operatorname{NB}(\mu_{i},k_{i})$. Then for any non-random weights
$w=(w_{1},\cdots,w_{n})^{\top}\in\mathbb{R}^{n}$ independent with
$\\{Y_{i}\\}_{i=1}^{n}$ and $t\geq 0$,

$\mathrm{P}\bigg{(}\Big{|}\sum_{i=1}^{n}w_{i}(Y_{i}-\mathrm{E}Y_{i})\Big{|}\geq
t\bigg{)}\leq
2\exp\left\\{-\frac{1}{4}\left(\frac{t^{2}}{2\sum_{i=1}^{n}w_{i}^{2}a^{2}(\mu_{i},k_{i})}\wedge\frac{t}{\displaystyle\max_{1\leq
i\leq n}|w_{i}|a(\mu_{i},k_{i})}\right)\right\\},$

where $q_{i}:=\frac{\mu_{i}}{k_{i}+\mu_{i}}\in(0,1)$ and
$a(\mu,k):=\left[\log\frac{1-(1-q)/\sqrt[k]{2}}{q}\right]^{-1}+\frac{\mu}{\log
2}.$

###### Proof.

We will use the sub-exponential norm. The moment generating function (MGF) for
$Y_{i}$ is

$\mathrm{E}e^{sY_{i}}=\left(\frac{1-q_{i}}{1-q_{i}e^{s}}\right)^{k_{i}}.$

Then by letting $\mathrm{E}\exp(|Y_{i}|/t)\leq 2$, we have

$2\geq\mathrm{E}\exp(|Y_{i}|/t)=\mathrm{E}\exp(Y_{i}/t)=\left(\frac{1-q_{i}}{1-q_{i}e^{1/t}}\right)^{k_{i}},$

which implies the sub-exponential norm for $Y_{i}$ is

$\|Y_{i}\|_{\psi_{1}}=\inf\\{t>0:\mathrm{E}\exp(|Y_{i}|/t)\leq
2\\}=\left[\log\frac{1-(1-q_{i})/\sqrt[k_{i}]{2}}{q_{i}}\right]^{-1}.$

Using the definition of $a_{i}$, from Proposition 4.2 in Zhang and Chen
(2021), we can immediately obtain the result in the Lemma. ∎

It should be noted that $a_{i}=a(\mu_{i},k_{i})$ naturally has a lower and
upper bound for any $i\in[n]$ since $\mu_{i}$ and $k_{i}$ are both bounded
between $e^{m_{s,n}}$ and $e^{M_{s,n}}.$

Note that $Y_{i}$ is an unbounded random variable; the next step is to find a
probabilistic bound for $M_{y,n}=\max_{i\in[n]}Y_{i}$. We will cite an
important lemma for this type of problem. We say a distribution $P_{\gamma}$
is strongly discrete log-concave with $\gamma>0$ if its density is strongly
midpoint log-convex with the same $\gamma>0$.

######  Lemma A.3.

[Concentration for strongly log-concave discrete distributions] Let
$P_{\gamma}$ be any strongly log-concave discrete distribution index by
$\gamma>0$ on $\mathbb{Z}^{n}$. Then for any function
$f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ that is L-Lipschitz with respect to
Euclidean norm, we have for $X\sim P_{\gamma}$,

$\mathrm{P}_{P_{\gamma}}\big{(}|f(X)-\mathrm{E}f(X)|\geq t\big{)}\leq
2\exp\left\\{-\frac{\gamma t^{2}}{4L^{2}}\right\\}$

for any $t>0$.

######  Lemma A.4.

The maximal of the response $M_{y,n}=\max_{i\in[n]}Y_{i}$ has the
concentration

$\mathrm{P}\left\\{M_{y,n}-\left(2\max_{i\in[n]}\left[a(\mu_{i},k_{i})-\frac{\mu_{i}}{\log
2}\right]\big{[}\log(2n)+\sqrt{2\log(2n)}\big{]}+\max_{i\in[n]}\mu_{i}\right)>t\right\\}\leq
e^{-\gamma t^{2}/4}$

for any $t>0$.

###### Proof.

For the upper bound of expectation, we first note that
$Y_{i}-\mathrm{E}Y_{i}\sim\operatorname{subE}(2\|Y_{i}\|_{\psi_{1}})$ with
$\|Y_{i}\|_{\psi_{1}}$ has calculated in Lemma A.2, then we have
$Y_{i}-\mathrm{E}Y_{i}\sim\operatorname{sub\Gamma}(4\|Y_{i}\|_{\psi_{1}}^{2},2\|Y_{i}\|_{\psi_{1}})$
by Example 5.3 in Zhang and Chen (2021), which further gives

$\displaystyle\mathrm{E}M_{y,n}$
$\displaystyle\leq\mathrm{E}\max_{i\in[n]}(Y_{i}-\mathrm{E}Y_{i})+\max_{i\in[n]}\mathrm{E}Y_{i}$
$\displaystyle\leq\big{(}2\cdot\max_{i\in[n]}4\|Y_{i}\|_{\psi_{1}}^{2}\cdot\log(2n)\big{)}^{\frac{1}{2}}+\max_{i\in[n]}2\|Y_{i}\|_{\psi_{1}}\cdot\log(2n)+\max_{i\in[n]}\mu_{i}$
$\displaystyle=2\max_{i\in[n]}\|Y_{i}\|_{\psi_{1}}\big{[}\log(2n)+\sqrt{2\log(2n)}\big{]}+\max_{i\in[n]}\mu_{i},$

where the second $\leq$ is by Corollary 7.3 in Zhang and Chen (2021) and the
bound in the lemma comes from the explicit expression in Lemma A.2.

By implementing Lemma A.3, the remaining we need to do is verifying
$\mathbb{Y}:=(Y_{1},\ldots,Y_{n})^{\top}\in\mathbb{Z}^{n}$ belongs to some
strongly log-concave discrete distribution $P_{\gamma}$ with the specifying
$\gamma>0$ after we take $f:(x_{1},\ldots,x_{n})\mapsto\max_{i\in[n]}x_{i}$
which is $1$-Lipschitz. By the definition, the derivative of log-density for
$y:=(y_{1},\ldots,y_{n})^{\top}$ is

$\psi^{\prime}(y_{i}):=\left.\frac{\partial\log p(y)}{\partial
y}\right|_{y_{i}}=\log\frac{\Gamma(k_{i}+y_{i})}{\Gamma(1+y_{i})}-y_{i}\log(k_{i}+\mu_{i}),$

then the Taylor expansion gives

$\displaystyle\psi(y)=\psi\left(\left[\frac{1}{2}x+\frac{1}{2}y\right]\right)+\frac{1}{2}\psi^{\prime}\left(\left[\frac{1}{2}x+\frac{1}{2}y\right]\right)(y-x)+\frac{1}{8}(y-x)^{2}\psi^{\prime\prime}\left(a_{1}\right),$
$\displaystyle\psi(x)=\psi\left(\left\lfloor\frac{1}{2}x+\frac{1}{2}y\right\rfloor\right)+\frac{1}{2}\psi^{\prime}\left(\left\lfloor\frac{1}{2}x+\frac{1}{2}y\right\rfloor\right)(x-y)+\frac{1}{8}(y-x)^{2}\psi^{\prime\prime}\left(a_{2}\right)$

where $a_{1}=t_{1}y+\left(1-t_{1}\right)(x+y)/2$,
$a_{2}=t_{2}y+\left(1-t_{1}\right)(x+y)/2$ with $t_{1},t_{2}\in[0,1]$. Define
the difference function

$\Delta(x,y):=\frac{x-y}{4}\left[\psi^{\prime}\left(\left[\frac{1}{2}x+\frac{1}{2}y\right]\right)-\psi^{\prime}\left(\left[\frac{1}{2}x+\frac{1}{2}y\right]\right)\right]+\frac{\psi^{\prime\prime}\left(a_{1}\right)+\psi^{\prime\prime}\left(a_{2}\right)}{16}(y-x)^{2},$

the Taylor expression above immediately implies

$\Delta(x,y)\geq|x-y|^{2}\left\\{\frac{\psi^{\prime\prime}\left(a_{1}\right)+\psi^{\prime\prime}\left(a_{2}\right)}{16}-\sup_{x\neq
y;x,y\in\mathbb{Z}^{n}}\frac{\left|\left[\psi^{\prime}(\lfloor(x+y)/2\rfloor)-\psi^{\prime}(\lceil(x+y)/2\rceil)\right]\right|}{4|x-y|}\right\\}.$

Let

$\displaystyle C_{\psi}$ $\displaystyle:=\sup_{x\neq
y;x,y\in\mathbb{Z}^{n}}\frac{\left|\left[\psi^{\prime}(\lfloor(x+y)/2\rfloor)-\psi^{\prime}(\lceil(x+y)/2\rceil)\right]\right|}{4|x-y|}$
$\displaystyle=\sup_{x\neq
y;x,y\in\mathbb{Z}^{n}}\left|\left[\log\frac{\Gamma(k_{i}+\lfloor(x+y)/2\rfloor)\Gamma(\lceil(x+y)/2\rceil+1)}{\Gamma(k_{i}+\lceil(x+y)/2\rceil)\Gamma(\lfloor(x+y)/2\rfloor+1)}-\frac{(\lfloor(x+y)/2\rfloor-\lceil(x+y)/2\rceil)}{\log^{-1}\left(k_{i}+\mu_{i}\right)}\right]\right|/4|x-y|,$

and it is not hard to see
$C_{\psi}\approx\frac{\left|\log\left(k_{i}+\mu_{i}\right)\right|}{4}$ or 0.
Besides,

$\displaystyle\psi^{\prime\prime}(y)$
$\displaystyle:=\left.\frac{\partial^{2}\log p(y)}{\partial
y^{2}}\right|_{y=y_{i}}=\frac{d}{dy_{i}}\log\frac{\Gamma\left(\theta+y_{i}\right)}{\Gamma\left(y_{i}+1\right)}$
$\displaystyle=\sum_{m=1}^{\infty}\left(\frac{1}{m+1}-\frac{1}{m+k_{i}+y_{i}}\right)-\sum_{m=1}^{\infty}\left(\frac{1}{m+1}-\frac{1}{m+y_{i}+1}\right)$
$\displaystyle=\sum_{m=1}^{\infty}\left(\frac{1}{m+y_{i}+1}-\frac{1}{m+k_{i}+y_{i}}\right)\geq\inf_{y_{i}\in\mathbb{Z}}\sum_{m=1}^{\infty}\left(\frac{1}{m+y_{i}+1}-\frac{1}{m+k_{i}+y_{i}}\right)=C_{\psi^{\prime\prime}}.$

Now, we have got

$\Delta(x,y)\geq|x-y|^{2}\left\\{\frac{\psi^{\prime\prime}\left(a_{1}\right)+\psi^{\prime\prime}\left(a_{2}\right)}{16}-C_{\psi}\right\\}\geq|x-y|^{2}\left(\frac{C_{\psi^{\prime\prime}}}{8}-C_{\psi}\right)$

which gives $\gamma=:\frac{C_{\psi^{\prime\prime}}}{8}-C_{\psi}>0$ from the
strong log-concave assumption for $\mathbb{Y}$, if
$C_{\psi}\approx\frac{|\log(k_{i}+\mu_{i})|}{4}$ is small. Hence, we can
conclude from Lemma A.3 and the upper bound of $\mathrm{E}M_{y,n}$

$\displaystyle\mathrm{P}\left\\{M_{y,n}-\left(2\max_{i\in[n]}\left[a(\mu_{i},k_{i})-\frac{\mu_{i}}{\log
2}\right]\big{[}\log(2n)+\sqrt{2\log(2n)}\big{]}+\max_{i\in[n]}\mu_{i}\right)>t\right\\}$
$\displaystyle\leq$ $\displaystyle\mathrm{P}(M_{y,n}-\mathrm{E}M_{y,n}>t)\leq
e^{-\gamma t^{2}/4}$

which is exact the result in the lemma. ∎

######  Remark A.1.

For $M_{y,n}$, it is distributed as sub-Gumbel, which is rarely studied by
research. Another way to deal with using the extreme value theory (EVT)
technique, we note that for any $t\in\mathbb{R}$

$\displaystyle\mathrm{P}(M_{y,n}-\mathrm{E}M_{y,n}>t)$
$\displaystyle=1-\prod_{i=1}^{n}\mathrm{P}\big{(}Y_{i}\leq
t+\mathrm{E}M_{y,n}\big{)}$
$\displaystyle=1-\prod_{i=1}^{n}\bigg{[}1-\exp\Big{\\{}-\frac{1}{4}\Big{(}\frac{(t+\mathrm{E}M_{y,n}-\mu_{i})^{2}}{2a_{i}^{2}}\wedge\frac{t+\mathrm{E}M_{y,n}-\mu_{i}}{a_{i}}\Big{)}\Big{\\}}\bigg{]}.$

If $Y_{i}$ is i.i.d, then in asymptotic sense,

$\displaystyle\mathrm{P}(M_{y,n}-\mathrm{E}M_{y,n}>t)$
$\displaystyle=1-\Big{[}1-\mathrm{P}\big{(}Y_{1}>t+\mathrm{E}M_{y,n}\big{)}\big{]}^{n}$
$\displaystyle\sim
1-\exp\left\\{-n\mathrm{P}\big{(}Y_{1}>t+\mathrm{E}M_{y,n}\big{)}\right\\}+o(1)$
$\displaystyle\sim
1-\exp\left\\{-n\exp\Big{(}-\frac{1}{4}\frac{(t+\mathrm{E}M_{y,n}-\mu_{1})^{2}}{2a_{1}^{2}}\wedge\frac{t+\mathrm{E}M_{y,n}-\mu_{1}}{a_{1}}\Big{)}\right\\}.$

Unfortunately, this technique cannot be used in the above lemma since: (i) we
need non-asymptotic version inequality instead of a vague expression with
$n\rightarrow\infty$; (ii) $\\{Y_{i}\\}$ is not an i.i.d series, and then EVT
theory will not be easily used in this particular setting. Hence we adopt a
discrete technique which has been use in Moriguchi et al. (2020) and fully
illustrated in Zhang and Jia (2022).

The stochastic Lipschitz conditions are established by using the properties of
$\partial_{1}\gamma(s,y)$ and $\partial_{2}\gamma(s,y)$. As we said before,
they are divided into two parts. The linear parts in them can be solved by the
concentration inequality for NB variables given in Lemma A.2, but the non-
linear part $\nu(s,y)$ needs some more advanced tools regarding the empirical
process. They are given as the following lemmas.

######  Lemma A.5 (The (3.12) in Sen (2018)).

Suppose $X_{1}(\omega),\cdots,X_{n}(\omega)\in\mathbb{R}$ are zero-mean
independent stochastic processes indexed by $\omega\in\Omega$. If there exist
$M_{0}$ and $S_{0}$ satisfying $|X_{i}(\omega)|\leq M_{0}$ and
$\sum_{i=1}^{n}\operatorname{var}\big{(}X_{i}(\omega)\big{)}\leq S_{0}^{2}$
for all $\omega\in\Omega$. Denote
$S_{n}=\sup_{\omega\in\Omega}\big{|}\sum_{i=1}^{n}X_{i}(\omega)\big{|}$, then
for any $t>0$,

$\mathrm{P}\big{(}S_{n}\geq 2\mathrm{E}S_{n}+S_{0}\sqrt{2t}+4M_{0}t\big{)}\leq
e^{-t}.$

A map $\phi:\mathbb{R}\rightarrow\mathbb{R}$ is called a contraction if
$|\phi(s)-\phi(t)|\leq|s-t|$ for all $s,t\in\mathbb{R}$. And in the following
lemmas $\varepsilon_{1},\cdots,\varepsilon_{n}$ are always i.i.d. Rademacher
variables.

######  Lemma A.6 (Theorem 2.2 in Chi (2010)).

Let $\mathcal{T}\subseteq V^{n}$ be a bounded set and $f_{1},\cdots,f_{n}$ be
functions $V\rightarrow\mathbb{R}$ such that $f_{i}$ is
$(M_{i},\ell_{\infty})$-Lipschitz with $f_{i}(0)=0$. For
$j=1,\cdots,k\in\mathbb{N}$, let
$T_{j}=\\{(t_{1j},\cdots,t_{nj}):(t_{1},\cdots,t_{n})\in\mathcal{T}\\}\subseteq\mathbb{R}^{n}$.
Then

$\mathrm{E}\sup_{t\in\mathcal{T}}\bigg{|}\sum_{i=1}^{n}\varepsilon_{i}f_{i}(t_{i})\bigg{|}\leq\beta_{k}\sum_{j=1}^{k}\mathrm{E}\sup_{s\in
T_{j}}\bigg{|}\sum_{i=1}^{n}\varepsilon_{i}M_{i}s_{i}\bigg{|},$

where $\beta_{k}$ is a universal constant that can be set no greater than
$3^{k}+3^{k-1}-2^{k}$.

######  Lemma A.7 (Theorem 4.12 in Ledoux and Talagrand (2013)).

Let $F:\mathbb{R}_{+}\rightarrow\mathbb{R}_{+}$ be convex and increasing. Let
further $\phi_{i}:\mathbb{R}\rightarrow\mathbb{R},i\leq n$ be contractions
such that $\phi_{i}(0)=0$. Then, for any bounded subset $\mathcal{T}$ in
$\mathbb{R}^{n}$,

$\mathrm{E}F\bigg{(}\frac{1}{2}\sup_{\mathcal{T}}\Big{|}\sum_{i=1}^{n}\varepsilon_{i}\phi_{i}(t_{i})\Big{|}\bigg{)}\leq\mathrm{E}F\bigg{(}\sup_{\mathcal{T}}\Big{|}\sum_{i=1}^{n}\varepsilon_{i}t_{i}\Big{|}\bigg{)}.$

######  Lemma A.8 (Lemma 5.2 in Massart (2000)).

Let $\mathcal{A}$ be some finite subset of $\mathbb{R}^{n}$, let
$R=\sup_{a\in\mathcal{A}}\left[\sum_{i=1}^{n}a_{i}^{2}\right]^{1/2}$, then

$\mathrm{E}\bigg{[}\sup_{a\in\mathcal{A}}\sum_{i=1}^{n}\varepsilon_{i}a_{i}\bigg{]}\leq
R\sqrt{2\log\big{(}\operatorname{card}(\mathcal{A})\big{)}}.$

With the assistance of these powerful tools, we can establish the stochastic
Lipschitz condition as follows, which is one of the most important points in
this article for establishing the oracle inequality of the $\ell_{2}$ distance
between the estimated value $\widehat{\theta}$ and the real value
$\theta^{*}$.

The proof of Corollary 3.1.

###### Proof.

Denote $c_{i}=Z_{i}\theta^{*}$. For $\theta\in\Theta$, denote
$t_{i}=Z_{i}(\theta-\theta^{*})=Z_{i}\theta-c_{i}$. We also define the map
$\bar{\pi}_{j}:(x_{1},\cdots,x_{p})^{\top}\mapsto(x_{1},\cdots,x_{j},0,\cdots,0)^{\top}$
and the function

$\varphi_{ij}(s)=\left\\{\begin{array}[]{ll}\displaystyle\frac{\gamma(c_{i}+\bar{\pi}_{j}s,Y_{i})-\gamma(c_{i}+\bar{\pi}_{j-1}s,Y_{i})}{s_{j}}-\partial_{j}\gamma(c_{i},Y_{i}),&\text{
if }s_{j}\neq 0;\\\
\partial_{j}\gamma(c_{i}+\bar{\pi}_{j-1}s,Y_{i})-\partial_{j}\gamma(c_{i},Y_{i}),&\text{
if }s_{j}=0.\end{array}\right.$

Thus, $\varphi_{ij}:\mathbb{R}^{2}\rightarrow\mathbb{R}$ is a real value
function for $i=1,\cdots,n,\,j=1,2$. Then it is easy to check that

$\gamma(Z_{i}\theta,Y_{i})-\gamma(Z_{i}\theta^{*},Y_{i})=\sum_{j=1}^{2}\big{(}\partial_{j}\gamma(c_{i},Y_{i})+\varphi_{ij}(t_{i})\big{)}t_{ij},$

and
$n\mathbb{P}_{n}\big{(}\gamma(\theta)-\gamma(\theta^{*})\big{)}=\sum_{i=1}^{n}\sum_{j=1}^{2}\big{(}\partial_{j}\gamma(c_{i},Y_{i})+\varphi_{ij}(t_{i})\big{)}X_{i}^{\top}\big{(}\theta^{(j)}-\theta^{*(j)}\big{)}$
in turn. It gives

$\displaystyle\sqrt{n}\mathbb{G}_{n}\big{(}\gamma(\theta)-\gamma(\theta^{*})\big{)}=$
$\displaystyle\sum_{i=1}^{n}\sum_{j=1}^{2}\big{(}\partial_{j}\gamma(c_{i},Y_{i})-\mathrm{E}\partial_{j}\gamma(c_{i},Y_{i})\big{)}X_{i}^{\top}\big{(}\theta^{(j)}-\theta^{*(j)}\big{)}$
$\displaystyle+\sum_{i=1}^{n}\sum_{j=1}^{2}\big{(}\varphi_{ij}(t_{i})-\mathrm{E}\varphi_{ij}(t_{i})\big{)}X_{i}^{\top}\big{(}\theta^{(j)}-\theta^{*(j)}\big{)}.$

First, we would like to give the explicit formula for $\varphi_{i1}$ and
obtain an upper bound as well as a Lipschitz parameter for $\varphi_{i2}$.
Denote $h_{i}(\cdot)=\gamma(\cdot,Y_{i})$, then

$\varphi_{ij}(s)=\int_{0}^{1}\big{(}\partial_{j}h_{i}(c_{i}+\bar{\pi}_{j-1}s+s_{j}ue_{j})-\partial_{j}h_{i}(c_{i})\big{)}\,du,$

where $e_{j}$ is the $j$-th basis vector of $\mathbb{R}^{2}$. Hence, for
$j=1$,

$\displaystyle\varphi_{i1}(s)$
$\displaystyle=-Y_{i}\int_{0}^{1}\left[\frac{e^{s_{2}}}{e^{s_{1}+u}+e^{s_{2}}}-\frac{e^{s_{2}}}{e^{s_{1}}+e^{s_{2}}}\right]\,du+\int_{0}^{1}\left[\frac{e^{s_{1}+s_{2}+u}}{e^{s_{1}+u}+e^{s_{2}}}-\frac{e^{s_{1}+s_{2}}}{e^{s_{1}}+e^{s_{2}}}\right]\,du$
$\displaystyle=\left[\log\frac{e^{s_{1}+1}+e^{s_{2}}}{e^{s_{1}}+e^{s_{2}}}-\frac{e^{s_{1}}}{e^{s_{1}}+e^{s_{2}}}\right]Y_{i}+C_{1}(s),$

in which $C_{1}(s)$ is a function only related to $s$ and free of $Y$ and the
index $i$. Using Lemma A.1, for $j=2$, write $F_{3}=F_{1}+M_{y,n}$,

$|\varphi_{i2}(s)|\leq\int_{0}^{1}\big{|}\partial_{2}h_{i}(c_{i}+\bar{\pi}_{1}s+s_{2}ue_{2})-\partial_{2}h_{i}(u)\big{|}\,du\leq
2F_{3},$

and

$\displaystyle\big{|}\varphi_{i2}(s)-\varphi_{i2}(t)\big{|}$
$\displaystyle\leq\int_{0}^{1}\big{|}\partial_{2}h_{i}(c_{i}+\bar{\pi}_{1}s+s_{2}ue_{2})-\partial_{2}h_{i}(c_{i}+\bar{\pi}_{1}t+t_{2}ue_{2})\big{|}\,du$
$\displaystyle\leq\int_{0}^{1}F_{2}\|\bar{\pi}_{1}(s-t)+(s_{2}-t_{2})ue_{2}\|_{\infty}\,du\leq
F_{2}\|s-t\|_{\infty}.$

This implies $\varphi_{i2}$ is $(F_{2},\ell_{\infty})$ Lipschitz. In
particular, letting $s=Z_{i}(\theta-\theta^{*})$ and $t=0$,

$\left|\varphi_{i2}\big{(}Z_{i}(\theta-\theta^{*})\big{)}\right|\leq\|Z_{i}(\theta-\theta^{*})\|_{\infty}\leq
F_{2}M_{x}D_{\Theta}.$

Hence, we obtain an upper bound for $\varphi_{i2}$ that

$\left|\varphi_{i2}\big{(}Z_{i}(\theta-\theta^{*})\big{)}\right|\leq
2F_{3}\vee F_{2}M_{x}D_{\Theta}:=M_{1}$ (A.1)

Now, for $k=1,\cdots,p$, define

$\xi_{ik}(\theta):=\big{(}\varphi_{i2}(t_{i})-\mathrm{E}\varphi_{i2}(t_{i})\big{)}X_{ik},\qquad
S_{k}=\sup_{\theta\in\Theta}\left|\sum_{i=1}^{n}\xi_{ik}(\theta)\right|.$

Then we can approach the final conclusion in the theorem by

$\displaystyle\sup_{\theta\in\Theta/\\{\theta^{*}\\}}\left|\frac{\sqrt{n}\mathbb{G}_{n}\big{(}\gamma(\theta)-\gamma(\theta^{*})\big{)}}{\|\theta-\theta^{*}\|_{1}}\right|$
$\displaystyle\leq\max_{1\leq k\leq
p}\left|\sum_{i=1}^{n}\big{(}\partial_{1}\gamma(c_{i},Y_{i})-\mathrm{E}\partial_{1}\gamma(c_{i},Y_{i})\big{)}X_{ik}\right|$
(A.2) $\displaystyle\qquad+\sup_{\theta\in\Theta/\\{\theta^{*}\\}}\max_{1\leq
k\leq
p}\left|\sum_{i=1}^{n}\big{(}\varphi_{i1}(t_{i})-\mathrm{E}\varphi_{i1}(t_{i})\big{)}X_{ik}\right|$
$\displaystyle\qquad+\max_{1\leq k\leq
p}\left|\sum_{i=1}^{n}\frac{e^{c_{i2}}}{e^{c_{i1}}+e^{c_{i2}}}\big{(}Y_{i}-\mathrm{E}Y_{i}\big{)}X_{ik}\right|$
$\displaystyle\qquad+\max_{1\leq k\leq
p}\left|\sum_{i=1}^{n}\big{(}\nu(c_{i},Y_{i})-\mathrm{E}\nu(c_{i},Y_{i})\big{)}X_{ik}\right|+\sup_{\theta\in\Theta/\\{\theta^{*}\\}}\max_{1\leq
k\leq p}S_{k}.$

We will tickle with (A.2) term by term.

_(i). the first three term in ( A.2):_

We will use concentration inequality to deal with these terms. For any $1\leq
k\leq p$ and $t\geq 0$, by Lemma A.2 and Cauchy-Schwartz inequality,

$\displaystyle\mathrm{P}\bigg{(}\Big{|}\sum_{i=1}^{n}\big{(}\partial_{1}\gamma(c_{i},Y_{i})-\mathrm{E}\partial_{1}\gamma(c_{i},Y_{i})\big{)}X_{ik}\Big{|}\geq
t\bigg{)}=\mathrm{P}\bigg{(}\Big{|}\frac{e^{c_{i2}}}{e^{c_{i1}}+e^{c_{i2}}}X_{ik}\big{(}Y_{i}-\mathrm{E}Y_{i}\big{)}\Big{|}\geq
t\bigg{)}$ $\displaystyle\leq
2\exp\left\\{-\frac{1}{4}\left(\frac{t^{2}}{2\sum_{i=1}^{n}(w_{i}^{(1)})^{2}X_{ik}^{2}a_{i}^{2}}\wedge\frac{t}{\displaystyle\max_{1\leq
i\leq n}|w_{i}^{(1)}X_{ik}|a_{i}}\right)\right\\}$ $\displaystyle\leq
2\exp\left\\{-\frac{1}{4}\left(\frac{t^{2}}{2\sqrt{\sum_{i=1}^{n}(w_{i}^{(1)})^{4}a_{i}^{4}}\displaystyle\max_{1\leq
k\leq
p}\sqrt{\begin{matrix}\sum_{i=1}^{n}X_{ik}^{4}\end{matrix}}}\wedge\frac{t}{M_{x}\displaystyle\max_{1\leq
i\leq n}|w_{i}^{(1)}|a_{i}}\right)\right\\},$

where $w_{i}^{(1)}=e^{c_{i2}}/(e^{c_{i1}}+e^{c_{i2}})$ and
$a_{i}=a(\mu_{i},k_{i})$ is defined in Lemma A.2, they are both determined and
free of $\theta$ and the index $k$. Hence,

$\displaystyle\mathrm{P}\bigg{(}\max_{1\leq k\leq
p}\Big{|}\sum_{i=1}^{n}\big{(}\partial_{1}\gamma(c_{i},Y_{i})-\mathrm{E}\partial_{1}\gamma(c_{i},Y_{i})\big{)}X_{ik}\Big{|}\geq
t\bigg{)}$ $\displaystyle\qquad\qquad\leq
2p\exp\left\\{-\frac{1}{4}\left(\frac{t^{2}}{2\sqrt{\sum_{i=1}^{n}(w_{i}^{(1)})^{4}a_{i}^{4}}\displaystyle\max_{1\leq
k\leq
p}\sqrt{\begin{matrix}\sum_{i=1}^{n}X_{ik}^{4}\end{matrix}}}\wedge\frac{t}{M_{x}\displaystyle\max_{1\leq
i\leq n}|w_{i}^{(1)}|a_{i}}\right)\right\\}.$

By letting the right side of the above display be $q_{1}\in(0,1)$, we can get

$\displaystyle\mathrm{P}\bigg{(}\max_{1\leq k\leq
p}\Big{|}\sum_{i=1}^{n}\big{(}\partial_{1}\gamma(c_{i},Y_{i})-\mathrm{E}\partial_{1}\gamma(c_{i},Y_{i})\big{)}X_{ik}\Big{|}$
$\displaystyle\quad\geq
2\sqrt{2\big{(}\sum_{i=1}^{n}(w_{i}^{(1)})^{4}a_{i}^{4}\big{)}^{1/2}\big{(}\max_{1\leq
k\leq p}\sum_{i=1}^{n}X_{ik}^{4}\big{)}^{1/2}\log(2p/q_{1})}\vee
4M_{x}\max_{1\leq i\leq n}|w_{i}^{(1)}|a_{i}\log(2p/q_{1})\Bigg{)}\leq q_{1}.$

Exactly the same, we can get for any $q_{3}\in(0,1)$, regarding to the third
term,

$\displaystyle\mathrm{P}\Bigg{(}\max_{1\leq k\leq
p}\Big{|}\sum_{i=1}^{n}\frac{e^{c_{i2}}}{e^{c_{i1}}+e^{c_{i2}}}\big{(}Y_{i}-\mathrm{E}Y_{i}\big{)}X_{ik}\Big{|}$
$\displaystyle\quad\geq
2\sqrt{2\big{(}\sum_{i=1}^{n}(w_{i}^{(1)})^{4}a_{i}^{4}\big{)}^{1/2}\big{(}\max_{1\leq
k\leq p}\sum_{i=1}^{n}X_{ik}^{4}\big{)}^{1/2}\log(2p/q_{3})}\vee
4M_{x}\max_{1\leq i\leq n}|w_{i}^{(1)}|a_{i}\log(2p/q_{3})\Bigg{)}\leq q_{3}.$

The situation is slightly different for the second term. Indeed,

$\displaystyle\mathrm{P}\left(\Big{|}\sum_{i=1}^{n}\big{(}\varphi_{i1}(t_{i})-\mathrm{E}\varphi_{i1}(t_{i})\big{)}X_{ik}\Big{|}\geq
t\right)$
$\displaystyle=\mathrm{P}\left(\left|\sum_{i=1}^{n}\left[\log\frac{e^{t_{i1}+1}+e^{t_{i2}}}{e^{t_{i1}}+e^{t_{i2}}}-\frac{e^{t_{i1}}}{e^{t_{i1}}+e^{t_{i2}}}\right]X_{ik}(Y_{i}-\mathrm{E}Y_{i})\right|\geq
t\right)$
$\displaystyle:=\mathrm{P}\left(\Big{|}\sum_{i=1}^{n}w_{i}^{(2)}(\theta)X_{ik}\big{(}Y_{i}-\mathrm{E}Y_{i}\big{)}\Big{|}\geq
t\right)$

Since $t_{i}$ is a function of $\theta$, so as the weights
$w_{i}^{(2)}(\theta)$, we cannot use the exact same method as previous.
However, since $\Theta$ is convex, we have
$\\{t_{i}\\}_{i=1}^{n}\subseteq\mathcal{S}$. Then it only needs to note that,

$\big{|}w_{i}^{(2)}(\theta)\big{|}=\left|\log\frac{e^{t_{i1}+1}+e^{t_{i2}}}{e^{t_{i1}}+e^{t_{i2}}}-\frac{e^{t_{i1}}}{e^{t_{i1}}+e^{t_{i2}}}\right|\leq\log\frac{e+e^{M_{s,n}-m_{s,n}}}{1+e^{m_{s,n}-M_{s,n}}}+\frac{1}{1+e^{m_{s,n}-M_{s,n}}}:=w^{(2)},$

which gives

$\displaystyle\mathrm{P}\Bigg{(}\max_{1\leq k\leq
p}\Big{|}\sum_{i=1}^{n}\big{(}\varphi_{i1}(t_{i})-\mathrm{E}\varphi_{i1}(t_{i})\big{)}X_{ik}\Big{|}$
$\displaystyle\geq
2\sqrt{2\sqrt{n}w^{(2)2}\big{(}\sum_{i=1}^{n}a_{i}^{4}\big{)}^{1/2}\big{(}\max_{1\leq
k\leq p}\sum_{i=1}^{n}X_{ik}^{4}\big{)}^{1/2}\log(2p/q_{2})}\vee
4M_{x}w^{(2)}\max_{1\leq i\leq n}|a_{i}|\log(2p/q_{2})\Bigg{)}\leq q_{2}.$

for any $\theta\in\Theta$ and $q_{2}\in(0,1)$.

_(ii). the forth term in ( A.2):_

From Lemma A.1, we know that $|\nu(c_{i},Y_{i})|\leq F_{1}$. Thus, simply by
Hoeffding inequality (see Corollary 2.1 (b) in Zhang and Chen (2021)), for any
$t\geq 0$ and $1\leq k\leq p$,

$\displaystyle\mathrm{P}\bigg{(}\Big{|}\sum_{i=1}^{n}\big{(}\nu(c_{i},Y_{i})-\mathrm{E}\nu(c_{i},Y_{i})\big{)}X_{ik}\Big{|}\geq
t\bigg{)}$ $\displaystyle\leq
2\exp\left\\{-\frac{t^{2}}{2F_{1}^{2}\sum_{i=1}^{n}X_{ik}^{2}}\right\\}\leq
2\exp\left\\{-\frac{t^{2}}{2F_{1}^{2}\displaystyle\max_{1\leq k\leq
p}\begin{matrix}\sum_{i=1}^{n}X_{ik}^{2}\end{matrix}}\right\\}.$

For arbitrary $q_{4}\in(0,1)$, let $t=F_{1}\sqrt{2\log(2p/q_{4})\max_{1\leq
k\leq p}\sum_{i=1}^{n}X_{ik}^{2}}$, we obtain

$\mathrm{P}\bigg{(}\max_{1\leq k\leq
p}\Big{|}\sum_{i=1}^{n}\big{(}\nu(c_{i},Y_{i})-\mathrm{E}\nu(c_{i},Y_{i})\big{)}X_{ik}\Big{|}\geq
F_{1}\sqrt{2\log(2p/q_{4})\max_{1\leq k\leq
p}\sum_{i=1}^{n}X_{ik}^{2}}\bigg{)}\leq q_{4}.$

_(iii). the last term in ( A.2):_

For any $i=1,\cdots,n$ and $k=1,\cdots,p$, by (A.1), $|\xi_{ik}(\theta)|\leq
2M_{1}M_{x}:=M_{0}$. And for any $\theta\in\Theta$, (A.1) also implies

$\sum_{i=1}^{n}\operatorname{var}\big{(}\xi_{ik}(\theta)\big{)}=\sum_{i=1}^{n}\mathrm{E}\big{(}\varphi_{i2}(t_{i})X_{ik}\big{)}^{2}\leq
A_{1}^{2}\sum_{i=1}^{n}X_{ik}^{2}\leq A_{1}^{2}\max_{1\leq k\leq
p}\sum_{i=1}^{n}X_{ik}^{2}:=S_{0}^{2}$

Therefore, from Lemma A.5, it follows that

$\mathrm{P}\left(S_{k}\geq 2\mathrm{E}S_{k}+S_{0}\sqrt{2t}+4M_{0}t\right)\leq
e^{-t}.$ (A.3)

So the last task is giving an upper bound for $\mathrm{E}S_{k}$. Note that
$\mathrm{E}\xi_{ik}(\theta)=0$, by symmetrization,

$\displaystyle\mathrm{E}S_{k}=\mathrm{E}\sup_{\theta\in\Theta}\bigg{|}\sum_{i=1}^{n}\big{(}\varphi_{i2}(t_{i})-\mathrm{E}\varphi_{i2}(t_{i})\big{)}X_{ik}\bigg{|}\leq
2\mathrm{E}\sup_{\theta\in\Theta}\bigg{|}\sum_{i=1}^{n}\varepsilon_{i}\varphi_{i2}(t_{i})X_{ik}\bigg{|}=2\mathrm{E}\sup_{t\in\mathcal{T}}\bigg{|}\sum_{i=1}^{n}\varepsilon_{i}\varphi_{i2}(t_{i})X_{ik}\bigg{|},$

where
$\mathcal{T}=\\{t_{i}=Z_{i}(\theta-\theta^{*}):\theta\in\Theta,i=1,\cdots,n\\}$,
and $\varepsilon_{1},\cdots,\varepsilon_{n}$ are i.i.d. Rademacher variables
independent of $Y_{1},\cdots,Y_{n}$. Here using the fact
$\varphi_{i2}(\cdot)X_{ik}$ is $(M_{x}F_{2},\ell_{\infty})$-Lipschitz and
Lemma A.6, A.7, and A.8,

$\displaystyle\mathrm{E}\sup_{t\in\mathcal{T}}\bigg{|}\sum_{i=1}^{n}\varepsilon_{i}\varphi_{i2}(t_{i})X_{ik}\bigg{|}$
$\displaystyle\leq
8M_{x}F_{2}\sum_{j=1}^{2}\mathrm{E}\sup_{t\in\mathcal{T}}\bigg{|}\varepsilon_{i}t_{ij}\bigg{|}=8M_{x}F_{2}\sum_{j=1}^{2}\mathrm{E}\sup_{\theta\in\Theta}\bigg{|}\sum_{i=1}^{n}\varepsilon_{i}X_{i}^{\top}\big{(}\theta^{(j)}-\theta^{*(j)}\big{)}\bigg{|}$
$\displaystyle\leq 16M_{x}F_{2}D_{\Theta}\mathrm{E}\max_{1\leq k\leq
p}\bigg{|}\sum_{i=1}^{n}\varepsilon_{i}X_{ik}\bigg{|}\leq 16\sqrt{2\log
p}M_{x}F_{2}D_{\Theta}\sqrt{\displaystyle\max_{1\leq k\leq
p}\sum_{i=1}^{n}X_{ik}^{2}}.$

Then by (A.3),

$\mathrm{P}\Bigg{(}S_{k}\geq 32\sqrt{2\log
p}M_{x}F_{2}D_{\Theta}\sqrt{\max_{1\leq k\leq
p}\sum_{i=1}^{n}X_{ik}^{2}}+M_{1}\sqrt{2t\max_{1\leq k\leq
p}\sum_{i=1}^{n}X_{ik}^{2}}+8M_{1}M_{x}t\Bigg{)}\leq e^{-t}.$

Note that the right side of the inequality is free of $\theta$, let
$t=\log(p/q_{5})$ in the above inequality, and use the same technique as
previous, we get the uniform bound for it. And the Theorem is proved by
letting $q_{2}=q_{3}=q_{1}$, $q_{4}=q_{2}$, $q_{5}=q_{3}$, and
$|w_{i}^{(1)}|\leq w^{(1)}$. ∎

The lower bound of the likelihood-based divergence

Recall the standard steps for establishing oracle inequality for lasso
estimator is (see Xiao et al. (2020) for example):

  * I.

To avoid the ill behavior of Hessian, propose the restricted eigenvalue
condition or other analogous conditions about the design matrix.

  * II.

Find the tuning parameter based on the high-probability event, i.e., the KKT
conditions.

  * III.

According to some restricted eigenvalue assumptions and tuning parameter
selection, derive the oracle inequalities via the definition of the Lasso
optimality and the minimizer under unknown expected risk function and some
basic inequalities. There are three sub-steps:

    * (i)

Under the KKT conditions, show that the error vector
$\widehat{\theta}-\theta^{*}$ is in some restricted set with structure
sparsity, and check that $\widehat{\theta}-\theta^{*}$ is in a big compact
set;

    * (ii)

Show that likelihood-based divergence of $\widehat{\theta}$ and $\theta^{*}$
can be lower bounded by some quadratic distance between $\widehat{\theta}$ and
$\theta^{*}$;

    * (iii)

By some elementary inequalities and (ii), show that
$\|\widehat{\theta}-\theta^{*}\|_{1}$ is in a smaller compact set with radius
of optimal rate (proportional to $\lambda$).

Under our approach, the KKT condition with a high probability is replaced by
the stochastic Lipschitz condition, while other steps should remain the same.
For most models belonging the canonical exponential family, the step III.(ii)
is quite trivial, see Lemma 1 in Abramovich and Grinshtein (2016) for example.
Nonetheless, it is worthy to note that the our loss function is not in the
canonical exponential family, so there is no extended discussion about the
lower bound of likelihood-based divergence of $\widehat{\theta}$ and
$\theta^{*}$ in our setting. We will use the following theorem to clarify this
thing.

######  Theorem A.1.

Suppose the condition is the same as that in Theorem 3.1. Denote the true
parameter for $Y_{i}$ is $\mu^{*}$ and $k^{*}$. If
$\\{Z_{i}\theta\\}_{i=1,\cdots,n,\theta\in\Theta}\subseteq\mathcal{S}\cap\\{s\in\mathbb{R}^{2}:2s_{1}+(1+s_{2}(1-k^{*})k^{*\mu^{*}})\mu^{*}\leq\frac{s_{1}+\mu^{*}}{2s_{2}^{2}}\\}$
and $\mu^{*}\geq 1$, then

$\mathrm{E}\gamma(Z_{i}\theta,Y_{i})-\mathrm{E}\gamma(Z_{i}\theta^{*},Y_{i})\geq
C_{\gamma}\|Z_{i}(\theta-\theta^{*})\|_{2}^{2},$

where $C_{\gamma}$ is a positive constant and its exact definition is in the
proof.

###### Proof.

For simplicity, we drop the index $i$. By the definition and the notation in
Theorem 3.1,

$\displaystyle\mathrm{E}\gamma(Z\theta,Y)-\mathrm{E}\gamma(Z\theta^{*},Y)$
$\displaystyle=D_{\operatorname{KL}}(s,c),$

where $D_{\operatorname{KL}}$ is the Kullback-Leibler divergence from the
$Y_{i}$’s density $f(y\,|\,Z\theta)$ to $f(y\,|\,Z\theta^{*})$, i.e.

$D_{\operatorname{KL}}(s,c):=\int
f(y\,|\,c)\log\frac{f(y\,|\,c)}{f(y\,|\,s)}\,dy.$

Due to the identification of the negative binomial distribution, we have
$D_{\operatorname{KL}}(s,c)\geq 0$ with equality if and only if $s=c$. Using
the Taylor theorem,

$\displaystyle D_{\operatorname{KL}}(s,c)$
$\displaystyle=D_{\operatorname{KL}}(c,c)+\left.\frac{\partial}{\partial
s}D_{\operatorname{KL}}(s,c)\right|_{s=c}+\frac{1}{2}(s-c)^{\top}\left[\frac{\partial^{2}}{\partial
s\partial s^{\top}}D_{\operatorname{KL}}(s,c)\right]_{s=c+\rho(s-c)}(s-c)$
$\displaystyle=\frac{1}{2}(s-c)^{\top}\left[\frac{\partial^{2}}{\partial
s\partial s^{\top}}D_{\operatorname{KL}}(s,c)\right]_{s=c+\rho(s-c)}(s-c)$
$\displaystyle\geq\frac{1}{2}\inf_{\rho\in[0,1]}\lambda_{min}\left[\frac{\partial^{2}}{\partial
s\partial
s^{\top}}D_{\operatorname{KL}}(s,c)\right]_{s=c+\rho(s-c)}\|s-c\|_{2}^{2}$

where $\rho\in[0,1]$ and $\lambda_{min}(\mathrm{M})$ is the smallest
eigenvalue of the matrix $\mathrm{M}$. So it is enough to show that
$\left[\frac{\partial^{2}}{\partial s\partial
s^{\top}}D_{\operatorname{KL}}(s,c)\right]_{s=c+\rho(s-c)}$ is strictly
positive define for any $\rho\in[0,1]$. First, calculate directly,

$\displaystyle\frac{\partial^{2}}{\partial s\partial
s^{\top}}D_{\operatorname{KL}}(s,c)$ $\displaystyle=\int
f(y\,|\,c)\left[\frac{\partial^{2}}{\partial s\partial
s^{\top}}\gamma(s,y)\right]\,dy$ $\displaystyle=\int
f(y\,|\,c)\left[\begin{array}[]{cc}\displaystyle\frac{e^{s_{1}+s_{2}}}{(e^{s_{1}}+e^{s_{2}})^{2}}(e^{s_{2}}+y)&\displaystyle\frac{e^{s_{1}+s_{2}}}{(e^{s_{1}}+e^{s_{2}})^{2}}(e^{s_{1}}-y)\\\
\displaystyle\frac{e^{s_{1}+s_{2}}}{(e^{s_{1}}+e^{s_{2}})^{2}}(e^{s_{1}}-y)&\displaystyle\partial_{2}\nu(s,y)+\frac{e^{s_{1}+s_{2}}}{(e^{s_{1}}+e^{s_{2}})^{2}}y\end{array}\right]\,dy$
$\displaystyle=:\left[\begin{array}[]{cc}a_{11}+b&a_{12}-b\\\
a_{21}-b&a_{22}+b\end{array}\right],$

where $a_{11}=\frac{e^{s_{1}+2s_{2}}}{(e^{s_{1}}+e^{s_{2}})^{2}},\qquad
a_{12}=\frac{e^{2s_{1}+s_{2}}}{(e^{s_{1}}+e^{s_{2}})^{2}},\qquad
b=\frac{e^{s_{1}+2s_{2}}}{(e^{s_{1}}+e^{s_{2}})^{2}}\mathrm{E}Y,$ and

$\displaystyle a_{22}=\mathrm{E}\partial_{2}v(s,Y)$
$\displaystyle=e^{s_{2}}\left[\psi(e^{s_{2}})+e^{s_{2}}\psi^{\prime}(e^{s_{2}})+\log(1+e^{s_{1}-s_{2}})-\frac{e^{s_{1}}}{e^{s_{1}}+e^{s_{2}}}-\Big{(}\frac{e^{s_{1}}}{e^{s_{1}}+e^{s_{2}}}\Big{)}^{2}\right]$
$\displaystyle\qquad\quad-e^{s_{2}}\left[\mathrm{E}\psi(Y+e^{s_{2}})+e^{s_{2}}\mathrm{E}\psi^{\prime}(Y+e^{s_{2}})\right].$

For a $2\times 2$ matrix $\mathrm{M}$, it is strictly positive define if and
only if $\operatorname{tr}(\mathrm{M})>0$ and
$\operatorname{det}(\mathrm{M})>0$. Denote $\mu=e^{s_{1}},k=e^{s_{2}}$, and
$\mu^{*}=e^{c_{1}},k^{*}=e^{c_{2}}$ are true parameters for $Y$. Then

$\displaystyle\operatorname{tr}\left[\frac{\partial^{2}}{\partial s\partial
s^{\top}}D_{\operatorname{KL}}(s,c)\right]$ $\displaystyle=\frac{\mu
k^{2}}{(\mu+k)^{2}}+2\frac{\mu
k^{2}}{(\mu+k)^{2}}\mu^{*}-k\left[\frac{\mu}{\mu+k}+\Big{(}\frac{\mu}{\mu+k}\Big{)}^{2}\right]$
(A.4)
$\displaystyle\qquad+k\left[\log(1+\mu/k)+\big{(}\psi(k)-\mathrm{E}\psi(Y+k)\big{)}+k\big{(}\psi^{\prime}(k)-\mathrm{E}\psi^{\prime}(Y+k)\big{)}\right]$
$\displaystyle=\frac{2(\mu^{*}-1)\mu
k^{2}}{(\mu+k)^{2}}+k\big{[}\log(1+\mu/k)+g_{1}(k)+kg_{2}(k)\big{]}$
$\displaystyle\geq k\big{[}\log(1+\mu/k)+g_{1}(k)+kg_{2}(k)\big{]}.$

Now, we are going to deal with $g_{1}(k)=\psi(k)-\mathrm{E}\psi(Y+k)$ and
$g_{2}(k)=\psi^{\prime}(k)-\mathrm{E}\psi^{\prime}(Y+k)$. For $\psi(x)$,

$0>\psi^{\prime\prime}(x)=-\frac{1}{x^{2}}-\int_{0}^{\infty}t^{2}\varphi(t)e^{-tx}\,dt\geq-\frac{1}{x^{2}}-\frac{2}{x^{3}}.$

Therefore, $\psi(\cdot)$ is concave. Using Jensen inequality and median value
theorem

$\displaystyle
g_{1}(k)=\psi(k)-\mathrm{E}\psi(Y+k)\geq\psi(k)-\psi\big{(}\mathrm{E}Y+k\big{)}\geq-\left(\frac{1}{k}+\frac{1}{k^{2}}\right)\mathrm{E}Y=-\mu^{*}\left(\frac{1}{k}+\frac{1}{k^{2}}\right).$

Similarly, for $g_{2}(k)$, by using the fact that
$\mathrm{E}(1/Y)=(1-k^{*})k^{*\mu^{*}}\mu^{*}$ and the assumption,

$\displaystyle g_{2}(k)$
$\displaystyle=\mathrm{E}\left[\psi^{\prime}(k)-\psi^{\prime}(Y+k)\right]\geq\mathrm{E}\left[Y\left(\frac{1}{(\xi(Y)+k)^{2}}+\frac{2}{(\xi(Y)+k)^{3}}\right)\right]$
$\displaystyle\geq\mathrm{E}\left[\frac{Y}{(Y+k)^{2}}\right]+2\mathrm{E}\left[\frac{Y}{(Y+k)^{3}}\right]\geq\left[\mathrm{E}\frac{(Y+k)^{2}}{Y}\right]^{-1}+2\left[\mathrm{E}\frac{(Y+k)^{3})}{Y}\right]^{-1}$
$\displaystyle=\frac{1}{2k+(1+k^{2}(1-k^{*})k^{*\mu^{*}})\mu^{*}}+\frac{2}{k^{*2}(k^{*}+\mu^{*})/\mu^{*2}+\mu^{*2}+3k\mu^{*}+3k^{2}+k^{3}(1-k^{*})k^{*\mu^{*}}\mu^{*}}$
$\displaystyle\geq(\mu+\mu^{*})\left(\frac{1}{2k^{2}}+\frac{1}{k^{3}}\right).$

where $\xi(Y)$ lies between $0$ and $Y$. The lower bounds for $g_{1}$ and
$g_{2}$, together with the fact that $\log(1+x)\geq x-x^{2}/2$ for $x\geq 0$,
we conclude that $\operatorname{tr}\left[\frac{\partial^{2}}{\partial
s\partial s^{\top}}D_{\operatorname{KL}}(s,c)\right]>0$. Similarly, we can
also prove $\operatorname{det}\left[\frac{\partial^{2}}{\partial s\partial
s^{\top}}D_{\operatorname{KL}}(s,c)\right]>0$, so the theorem holds. ∎

The proof of Theorem 3.1

###### Proof.

The proof follows the idea in Bickel et al. (2009). First, by the definition
of $\widehat{\theta}$,

$\displaystyle P\big{(}\gamma(\widehat{\theta})-\gamma(\theta^{*})\big{)}$
$\displaystyle\leq
P\big{(}\gamma(\widehat{\theta})-\gamma(\theta^{*})\big{)}+\big{(}\mathbb{P}_{n}\gamma(\theta^{*})+\lambda\|\theta^{*}\|_{\omega,1}\big{)}-\big{(}\mathbb{P}_{n}\gamma(\widehat{\theta})+\lambda\|\widehat{\theta}\|_{\omega,1}\big{)}$
$\displaystyle\leq\frac{1}{\sqrt{n}}\mathbb{G}_{n}\big{(}\gamma(\theta^{*})-\gamma(\widehat{\theta})\big{)}+\lambda\big{(}\|\theta^{*}\|_{\omega,1}-\|\widehat{\theta}\|_{\omega,1}\big{)}.$

From Theorem A.1, we also have

$P\big{(}\gamma(\widehat{\theta})-\gamma(\theta^{*})\big{)}\geq\frac{C_{\gamma}}{n}\sum_{i=1}^{n}\|Z_{i}(\widehat{\theta}-\theta^{*})\|_{2}^{2}=\frac{C_{\gamma}}{n}\sum_{j=1}^{2}\|\mathrm{X}(\widehat{\theta}^{(j)}-\theta^{*(j)})\|_{2}^{2}.$

Then by Theorem 3.1 and the definition of $\lambda$,

$\displaystyle
C_{\gamma}\sum_{j=1}^{2}\|\mathrm{X}(\widehat{\theta}^{(j)}-\theta^{*(j)})\|_{2}^{2}$
$\displaystyle\leq{\sqrt{n}}\mathbb{G}_{n}\big{(}\gamma(\theta^{*})-\gamma(\widehat{\theta})\big{)}+n\lambda\big{(}\|\theta^{*}\|_{\omega,1}-\|\widehat{\theta}\|_{\omega,1}\big{)}$
$\displaystyle\leq
M_{q}\|\theta^{*}-\widehat{\theta}\|_{1}+(1+1/a)M_{q}\big{(}\|\theta^{*}\|_{\omega,1}-\|\widehat{\theta}\|_{\omega,1}\big{)}$
$\displaystyle=M_{q}\sum_{j=1}^{2}\bigg{[}\|\widehat{\theta}^{(j)}-\theta^{*(j)}\|_{1}+(1+1/a)\omega_{j}\big{(}\|\theta^{*(j)}\|_{1}-\|\widehat{\theta}^{(j)}\|_{1}\big{)}\bigg{]}$

holds with probability at least $1-q$, where $a=(K-1)/2$. Now let
$J_{1},J_{2}\subseteq\\{1,\cdots,p\\}$ be any sets with
$J_{j}\supseteq\operatorname{spt}\big{(}\theta^{*(j)}\big{)}$. It is easy to
check

$\displaystyle\|\widehat{\theta}^{(j)}-\theta^{*(j)}\|_{1}$
$\displaystyle+(1+1/a)\omega_{j}\big{(}\|\theta^{*(j)}\|_{1}-\|\widehat{\theta}^{(j)}\|_{1}\big{)}$
$\displaystyle=\|\widehat{\theta}_{J_{j}}^{(j)}-\theta^{*(j)}\|_{1}+\|\widehat{\theta}_{J_{j}^{c}}^{(j)}\|_{1}+(1+1/a)\omega_{j}\big{(}\|\theta^{*(j)}\|_{1}-\|\widehat{\theta}_{J_{j}}^{(j)}\|_{1}-\|\widehat{\theta}_{J_{j}^{c}}^{(j)}\|_{1}\big{)}$
$\displaystyle\leq(K/a)\|\widehat{\theta}_{J_{j}}^{(j)}-\theta^{*(j)}\|_{1}-(1/a)\|\widehat{\theta}_{J_{j}^{c}}^{(j)}\|_{1}.$

by the fact $\omega_{j}\in[0,1]$. It gives that with probability at least
$1-q$,

$\sum_{j=1}^{2}\|\mathrm{X}(\widehat{\theta}^{(j)}-\theta^{*(j)})\|_{2}^{2}\leq\frac{M_{q}}{aC_{\gamma}}\sum_{j=1}^{2}\big{(}K\|\widehat{\theta}_{J_{j}}^{(j)}-\theta^{*(j)}\|_{1}-\|\widehat{\theta}_{J_{j}^{c}}^{(j)}\|_{1}\big{)}.$
(A.5)

Let $A_{1},A_{2}\subseteq\\{1,\cdots,p\\}$ satisfying
$\operatorname{spt}\big{(}\theta^{*(j)}\big{)}\subseteq A_{j}$ and
$\operatorname{card}(A_{j})=p_{1}$, and we also let $B_{j}$ be the union of
$A_{j}$ and the indices of $p_{1}$ largest $\widehat{\theta}^{(j)}$. Then
$A_{j}$ and $B_{j}$ also guarantee (A.5). And from Lemma 3.1, they also give

$\|\widehat{\theta}^{(j)}_{B_{j}^{c}}\|_{2}^{2}\leq
p_{1}^{-1}\|\widehat{\theta}^{(j)}_{A_{j}^{c}}\|_{1}^{2}.$

Besides, from the definition of $A_{j}$ and $B_{j}$, we know that
$\|\widehat{\theta}^{(j)}_{A_{j}^{c}}\|_{1}\geq\|\widehat{\theta}^{(j)}_{B_{j}^{c}}\|_{1}$
and
$\|\widehat{\theta}^{(j)}_{A_{j}}-\theta^{*(j)}\|_{1}\leq\|\widehat{\theta}^{(j)}_{B_{j}}-\theta^{*(j)}\|_{1}$.

Unlike the single lasso question, here we need to define
$I:=\\{j=1,2:K\|\widehat{\theta}_{A_{j}}^{(j)}-\theta^{*(j)}\|_{1}\geq\|\widehat{\theta}_{A_{j}^{c}}^{(j)}\|_{1}\\}$,
and consider $j\in I$ and $j\notin I$ separately. Obviously, $I\neq\emptyset$,
or (A.5) cannot be holden. For $j\in I$, we have

$K\|\widehat{\theta}^{(j)}_{B_{j}}-\theta^{*(j)}\|_{1}-\|\widehat{\theta}^{(j)}_{B_{j}^{c}}\|_{1}\geq
K\|\widehat{\theta}^{(j)}_{A_{j}}-\theta^{*(j)}\|_{1}-\|\widehat{\theta}^{(j)}_{A_{j}^{c}}\|_{1}\geq
0.$

Then by the restricted eigenvalue condition,

$n\kappa^{2}\|\widehat{\theta}_{J_{j}}^{(j)}-\theta^{*(j)}\|_{2}^{2}\leq\|\mathrm{X}(\widehat{\theta}^{(j)}-\theta^{*(j)})\|_{2}^{2}$

holds for $J_{j}=A_{j}$ or $J_{j}=B_{j}$. Note that from (A.5),

$\sum_{j\in
I}\|\mathrm{X}(\widehat{\theta}^{(j)}-\theta^{*(j)})\|_{2}^{2}\leq\frac{M_{q}}{aC_{\gamma}}\sum_{j\in
I}\big{(}\|\widehat{\theta}_{A_{j}}^{(j)}-\theta^{*(j)}\|_{1}-\|\widehat{\theta}_{A_{j}^{c}}^{(j)}\|_{1}\big{)}\leq\frac{M_{q}}{aC_{\gamma}}\sum_{j\in
I}\big{(}\|\widehat{\theta}_{B_{j}}^{(j)}-\theta^{*(j)}\|_{1}-\|\widehat{\theta}_{B_{j}^{c}}^{(j)}\|_{1}\big{)},$

then by Cauchy-Schwartz inequality,

$\displaystyle n\kappa^{2}\sum_{j\in
I}\|\widehat{\theta}_{A_{j}}^{(j)}-\theta^{*(j)}\|_{2}^{2}$
$\displaystyle\leq\|\mathrm{X}(\widehat{\theta}_{A_{j}}^{(j)}-\theta^{*(j)})\|_{2}^{2}\leq\frac{M_{q}K}{aC_{\gamma}}\sum_{j\in
I}\|\widehat{\theta}_{A_{j}}^{(j)}-\theta^{*(j)}\|_{1}$
$\displaystyle\leq\frac{M_{q}K\sqrt{p_{1}}}{aC_{\gamma}}\sum_{j\in
I}\|\widehat{\theta}_{A_{j}}^{(j)}-\theta^{*(j)}\|_{2}\leq\frac{M_{q}K\sqrt{2p_{1}}}{aC_{\gamma}}\left[\sum_{j\in
I}\|\widehat{\theta}_{A_{j}}^{(j)}-\theta^{*(j)}\|_{2}^{2}\right]^{1/2}.$

It gives

$\sum_{j\in
I}\|\widehat{\theta}_{A_{j}}^{(j)}-\theta^{*(j)}\|_{2}^{2}\leq\frac{2p_{1}M_{q}^{2}K^{2}}{a^{2}\kappa^{4}n^{2}C_{\gamma}^{2}},\qquad\sum_{j\in
I}\|\widehat{\theta}_{B_{j}}^{(j)}-\theta^{*(j)}\|_{2}^{2}\leq\frac{4p_{1}M_{q}^{2}K^{2}}{a^{2}\kappa^{4}n^{2}C_{\gamma}^{2}},$

where we use that fact $\operatorname{card}(B_{j})=2p_{1}$. Furthermore, since

$\displaystyle\|\widehat{\theta}^{(j)}_{B_{j}^{c}}\|_{2}^{2}\leq\sum_{j\in
I}p_{1}^{-1}{\|\widehat{\theta}^{(j)}_{A_{j}^{c}}\|_{1}^{2}}\leq\frac{K^{2}}{p_{1}}\sum_{j\in
I}{\|\widehat{\theta}_{A_{j}}^{(j)}-\theta^{*(j)}\|_{1}^{2}}\leq
K^{2}\sum_{j\in I}{\|\widehat{\theta}_{A_{j}}^{(j)}-\theta^{*(j)}\|_{2}^{2}},$

we can conclude that

$\displaystyle\sum_{j\in I}\|\widehat{\theta}^{(j)}-\theta^{*(j)}\|_{2}^{2}$
$\displaystyle=\sum_{j\in
I}\big{(}\|\widehat{\theta}_{B_{j}}^{(j)}-\theta^{*(j)}\|_{2}^{2}+\|\widehat{\theta}^{(j)}_{B_{j}^{c}}\|_{2}^{2}\big{)}$
(A.6) $\displaystyle\leq\sum_{j\in
I}\big{(}\|\widehat{\theta}_{B_{j}}^{(j)}-\theta^{*(j)}\|_{2}^{2}+K^{2}\|\widehat{\theta}_{A_{j}}^{(j)}-\theta^{*(j)}\|_{2}^{2}\big{)}=\frac{2p_{1}M_{q}^{2}(2+K^{2})K^{2}}{a^{2}\kappa^{4}n^{2}C_{\gamma}^{2}}.$

Now we will tickle the situation that $j\notin I$. For $j\notin I$,
$K\|\widehat{\theta}_{A_{j}}^{(j)}-\theta^{*(j)}\|_{1}<\|\widehat{\theta}_{A_{j}^{c}}^{(j)}\|_{1}$.
Again from (A.5), we have

$\sum_{j\notin
I}\|\mathrm{X}(\widehat{\theta}^{(j)}-\theta^{*(j)})\|_{2}^{2}\leq\frac{M_{q}K}{aC_{\gamma}}\sum_{j\in
I}\|\widehat{\theta}_{A_{j}}^{(j)}-\theta^{*(j)}\|_{1}$

and

$0\leq\sum_{j\notin
I}\big{(}\|\widehat{\theta}_{A_{j}^{c}}^{(j)}\|_{1}-K\|\widehat{\theta}_{A_{j}}^{(j)}-\theta^{*(j)}\|_{1}\big{)}\leq
K\sum_{j\in I}\|\widehat{\theta}_{A_{j}}^{(j)}-\theta^{*(j)}\|_{1}.$

Indeed, if the two inequalities above has the opposite direction, then for the
first one, one can find that

$\sum_{j\in
I}\|\mathrm{X}(\widehat{\theta}^{(j)}-\theta^{*(j)})\|_{2}^{2}\leq\frac{M_{q}}{aC_{\gamma}}\left[\sum_{j\notin
I}\big{(}K\|\widehat{\theta}_{A_{j}}^{(j)}-\theta^{*(j)}\|_{1}-\|\widehat{\theta}_{A_{j}^{c}}^{(j)}\|_{1}\big{)}-\sum_{j\in
I}\|\widehat{\theta}_{A_{j}^{c}}^{(j)}\|_{1}\right]<0,$

and

$\sum_{j=1}^{2}\|\mathrm{X}(\widehat{\theta}^{(j)}-\theta^{*(j)})\|_{2}^{2}\leq-\frac{M_{q}}{aC_{\gamma}}\sum_{j\in
I}\|\widehat{\theta}_{A_{j}^{c}}^{(j)}\|_{1}<0.$

Once again, by Cauchy-Schwartz inequality,

$\sum_{j\in
I}\|\widehat{\theta}_{A_{j}}^{(j)}-\theta^{*(j)}\|_{1}\leq\sqrt{p_{1}}\sum_{j\in
I}\|\widehat{\theta}_{A_{j}}^{(j)}-\theta^{*(j)}\|_{2}\leq\sqrt{2p_{1}}\left[\sum_{j\in
I}\|\widehat{\theta}_{A_{j}}^{(j)}-\theta^{*(j)}\|_{2}^{2}\right]^{1/2}\leq\frac{2p_{1}M_{q}K}{a\kappa^{2}nC_{\gamma}}.$

Denote
$\Delta_{j}:=\|\widehat{\theta}_{A_{j}^{c}}^{(j)}\|_{1}-K\|\widehat{\theta}_{A_{j}}^{(j)}-\theta^{*(j)}\|_{1}$.
Then for $j\notin J$, $\Delta_{j}>0$, and

$\sum_{j\notin I}\Delta_{j}\leq K\sum_{j\in
I}\|\widehat{\theta}_{A_{j}}^{(j)}-\theta^{*(j)}\|_{1}\leq\frac{2p_{1}M_{q}K^{2}}{a\kappa^{2}nC_{\gamma}}.$

For any $j\notin I$, define

$\widetilde{\theta}^{(j)}=\widehat{\theta}^{(j)}+\frac{\Delta_{j}}{p_{1}K}\sum_{k\in
A_{j}}\operatorname{sgn}\big{(}\widehat{\theta}_{k}^{(j)}-\theta_{k}^{*(j)}\big{)}e_{k}.$

Then for $k\in A_{j}$,

$|\widetilde{\theta}_{k}^{(j)}-\theta_{k}^{*(j)}|=|\widehat{\theta}_{k}^{(j)}-\theta_{k}^{*(j)}|+\frac{\Delta_{j}}{p_{1}K},$

while for $k\notin I$,
$\widetilde{\theta}^{(j)}_{k}=\widehat{\theta}^{(j)}_{k}$. Therefore,

$K\|\widetilde{\theta}_{A_{j}}^{(j)}-\theta^{*(j)}\|_{1}=K\left[\|\widehat{\theta}_{A_{j}}^{(j)}-\theta^{*(j)}\|_{1}+\sum_{k\in
A_{j}}\frac{\Delta_{j}}{p_{1}K}\right]=\|\widehat{\theta}_{A_{j}^{c}}^{(j)}\|_{1}=\|\widetilde{\theta}_{A_{j}^{c}}^{(j)}\|_{1},$

and consequently $\|\widetilde{\theta}_{B_{j}^{c}}^{(j)}\|_{1}\leq
K\|\widetilde{\theta}_{B_{j}}^{(j)}-\theta^{*(j)}\|_{1}$. Once again by the
restricted eigenvalue condition,

$\|\mathrm{X}(\widetilde{\theta}^{(j)}-\theta^{*(j)})\|_{2}^{2}\geq
n\kappa^{2}\|\widetilde{\theta}_{B_{j}}^{(j)}-\theta^{*(j)}\|_{2}^{2}\geq
n\kappa^{2}\|\widetilde{\theta}_{A_{j}}^{(j)}-\theta^{*(j)}\|_{2}^{2}.$ (A.7)

On the other hand, note that for any $s,t\in\mathbb{R}^{m}$ inequality
$\|s+t\|_{2}^{2}\leq 2(\|s\|_{2}^{2}+\|t\|_{2}^{2})$ and
$\|s\|_{2}\leq\|s\|_{1}\leq\sqrt{m}\|s\|_{2}$ hold, we conclude

$\displaystyle\sum_{j\notin
I}\|\mathrm{X}(\widetilde{\theta}^{(j)}-\theta^{*(j)})\|_{2}^{2}$
$\displaystyle\leq 2\sum_{j\notin
I}\Big{(}\|\mathrm{X}(\widehat{\theta}^{(j)}-\theta^{*(j)})\|_{2}^{2}+\|\mathrm{X}(\widehat{\theta}^{(j)}-\widetilde{\theta}^{(j)})\|_{2}^{2}\Big{)}$
(A.8) $\displaystyle\leq\frac{2M_{q}K}{aC_{\gamma}}\sum_{j\in
I}\|\widehat{\theta}_{A_{j}}^{(j)}-\theta^{*(j)}\|_{1}+2\sum_{j\notin
I}\|\mathrm{X}(\widehat{\theta}^{(j)}-\widetilde{\theta}^{(j)})\|_{2}^{2}$
$\displaystyle\leq\frac{4p_{1}M_{q}^{2}K^{2}}{na^{2}\kappa^{2}C_{\gamma}^{2}}+2\sum_{j\notin
I}\|\mathrm{X}(\widehat{\theta}^{(j)}-\widetilde{\theta}^{(j)})\|_{2}^{2}.$

Next, we will use the definition of the $p_{1}$-restricted isometry constant
$\sigma_{\mathrm{X},l}^{2}$. Since
$\operatorname{spt}\big{(}\widetilde{\theta}^{(j)}-\widehat{\theta}^{(j)}\big{)}\leq\operatorname{card}(A_{j})=p_{1}$,
then

$\displaystyle\sum_{j\notin
I}\|\mathrm{X}(\widehat{\theta}^{(j)}-\widetilde{\theta}^{(j)})\|_{2}^{2}$
$\displaystyle\leq\sigma_{\mathrm{X},p_{1}}^{2}\sum_{j\notin
I}\|\widehat{\theta}^{(j)}-\widetilde{\theta}^{(j)}\|_{2}^{2}$
$\displaystyle=\sigma_{\mathrm{X},p_{1}}^{2}\sum_{j\notin I}\sum_{k\in
A_{j}}\Big{(}\frac{\Delta_{j}}{p_{1}K}\Big{)}^{2}=\frac{\sigma_{\mathrm{X},p_{1}}^{2}}{p_{1}K^{2}}\sum_{j\notin
I}\Delta_{j}^{2}$
$\displaystyle\leq\frac{\sigma_{\mathrm{X},p_{1}}^{2}}{p_{1}K^{2}}\Big{(}\sum_{j\notin
I}\Delta_{j}\Big{)}^{2}\leq\frac{4p_{1}\sigma_{\mathrm{X},p_{1}}^{2}K^{2}}{a^{2}\kappa^{4}n^{2}C_{\gamma}^{2}}.$

The above inequality together with (A.7) and (A.8) gives

$\sum_{j\notin
I}\|\widetilde{\theta}_{A_{j}}^{(j)}-\theta^{*(j)}\|_{2}^{2}\leq\sum_{j\notin
I}\|\widetilde{\theta}_{B_{j}}^{(j)}-\theta^{*(j)}\|_{2}^{2}\leq\frac{4p_{1}(n\kappa^{2}+2\sigma_{\mathrm{X},p_{1}}^{2})M_{q}^{2}K^{2}}{a^{2}C_{\gamma}^{2}n^{3}\kappa^{6}}.$

Finally, since

$\|\widetilde{\theta}_{B_{j}^{c}}^{(j)}\|_{2}^{2}\leq\|\widetilde{\theta}_{B_{j}^{c}}^{(j)}\|_{1}^{2}\leq
K^{2}\|\widetilde{\theta}_{B_{j}}^{(j)}-\theta^{*(j)}\|_{1}^{2}\leq
2p_{1}K^{2}\|\widetilde{\theta}_{B_{j}}^{(j)}-\theta^{*(j)}\|_{2}^{2},$

we obtain that

$\displaystyle\sum_{j\notin
I}\|\widehat{\theta}^{(j)}-\theta^{*(j)}\|_{2}^{2}$
$\displaystyle\leq\sum_{j\notin
I}\|\widetilde{\theta}^{(j)}-\theta^{*(j)}\|_{2}^{2}=\sum_{j\notin
I}\big{(}\|\widetilde{\theta}_{B_{j}}^{(j)}-\theta^{*(j)}\|_{2}^{2}+\|\widetilde{\theta}_{B_{j}^{c}}^{(j)}\|_{2}^{2}\big{)}$
(A.9) $\displaystyle\leq(1+2p_{1}K)\sum_{j\notin
I}\|\widetilde{\theta}_{B_{j}}^{(j)}-\theta^{*(j)}\|_{2}^{2}\leq\frac{4p_{1}(1+2p_{1}K)(n\kappa^{2}+2\sigma_{\mathrm{X},p_{1}}^{2})M_{q}^{2}K^{2}}{a^{2}C_{\gamma}^{2}n^{3}\kappa^{6}}.$

Combining (A.6) and (A.9), it is easy to see the remaining ∎

## References

  * Abramovich and Grinshtein (2016) Abramovich, F. and V. Grinshtein (2016).  Model selection and minimax estimation in generalized linear models.  IEEE Transactions on Information Theory 62(6), 3721–3730.
  * Adamczak (2008) Adamczak, R. (2008).  A tail inequality for suprema of unbounded empirical processes with applications to markov chains.  Electronic Journal of Probability 13, 1000–1034.
  * Alzer (1997) Alzer, H. (1997).  On some inequalities for the gamma and psi functions.  Mathematics of computation 66(217), 373–389.
  * Bateman (1953) Bateman, H. (1953).  Higher transcendental functions [volumes i-iii], Volume 1.  McGraw-Hill Book Company.
  * Bickel et al. (2009) Bickel, P. J., Y. Ritov, A. B. Tsybakov, et al. (2009).  Simultaneous analysis of lasso and dantzig selector.  The Annals of statistics 37(4), 1705–1732.
  * Candes et al. (2007) Candes, E., T. Tao, et al. (2007).  The dantzig selector: Statistical estimation when p is much larger than n.  The annals of Statistics 35(6), 2313–2351.
  * Chi (2010) Chi, Z. (2010).  Stochastic lipschitz continuity for high dimensional lasso with multiple linear covariate structures or hidden linear covariates.  arXiv preprint arXiv:1011.1384.
  * Cui et al. (2021) Cui, C., J. Jia, Y. Xiao, and H. Zhang (2021).  Directional fdr control for sub-gaussian sparse glms.  arXiv preprint arXiv:2105.00393.
  * Dai et al. (2013) Dai, H., Y. Bao, and M. Bao (2013).  Maximum likelihood estimate for the dispersion parameter of the negative binomial distribution.  Statistics & Probability Letters 83(1), 21–27.
  * Hilbe (2011) Hilbe, J. M. (2011).  Negative binomial regression.  Cambridge University Press.
  * Huang et al. (2021) Huang, H., H. Zhang, and B. Li (2021).  Weighted lasso estimates for sparse logistic regression: non-asymptotic properties with measurement errors.  Acta Mathematica Scientia 41(1), 207–230.
  * Ledoux and Talagrand (2013) Ledoux, M. and M. Talagrand (2013).  Probability in Banach Spaces: isoperimetry and processes.  Springer Science & Business Media.
  * Massart (2000) Massart, P. (2000).  Some applications of concentration inequalities to statistics.  In Annales de la Faculté des sciences de Toulouse: Mathématiques, Volume 9, pp. 245–303.
  * Moriguchi et al. (2020) Moriguchi, S., K. Murota, A. Tamura, and F. Tardella (2020).  Discrete midpoint convexity.  Mathematics of Operations Research 45(1), 99–128.
  * Nguelifack and Kemajou-Brown (2019) Nguelifack, B. M. and I. Kemajou-Brown (2019).  Robust rank-based variable selection in double generalized linear models with diverging number of parameters under adaptive lasso.  Journal of Statistical Computation and Simulation 89(11), 2051–2072.
  * Qiu et al. (2018) Qiu, Y., S. X. Chen, and D. Nettleton (2018).  Detecting rare and faint signals via thresholding maximum likelihood estimators.  The Annals of Statistics 46(2), 895–923.
  * Riphahn et al. (2003) Riphahn, R. T., A. Wambach, and A. Million (2003).  Incentive effects in the demand for health care: a bivariate panel count data estimation.  Journal of applied econometrics 18(4), 387–405.
  * Sen (2018) Sen, B. (2018).  A gentle introduction to empirical process theory and applications.
  * Shi et al. (2019) Shi, C., R. Song, Z. Chen, and R. Li (2019).  Linear hypothesis testing for high dimensional generalized linear models.  Annals of statistics 47(5), 2671.
  * Tian et al. (2020) Tian, S., D. Wang, and S. Cui (2020).  A seasonal geometric inar process based on negative binomial thinning operator.  Statistical Papers 61(6), 2561–2581.
  * Tibshirani (1996) Tibshirani, R. (1996).  Regression shrinkage and selection via the lasso.  Journal of the Royal Statal Society, Series B 58(1).
  * Wang et al. (2016) Wang, Z., S. Ma, M. Zappitelli, C. Parikh, C.-Y. Wang, and P. Devarajan (2016).  Penalized count data regression with application to hospital stay after pediatric cardiac surgery.  Statistical methods in medical research 25(6), 2685–2703.
  * Weißbach and Radloff (2020) Weißbach, R. and L. Radloff (2020).  Consistency for the negative binomial regression with fixed covariate.  Metrika 83(5), 627–641.
  * Xiao et al. (2020) Xiao, Y., T. Yan, H. Zhang, and Y. Zhang (2020).  Oracle inequalities for weighted group lasso in high-dimensional misspecified cox models.  Journal of Inequalities and Applications 2020(1), 1–33.
  * Xie and Lederer (2021) Xie, F. and J. Lederer (2021).  Aggregating knockoffs for false discovery rate control with an application to gut microbiome data.  Entropy 23(2), 230.
  * Xie and Xiao (2020) Xie, F. and Z. Xiao (2020).  Consistency of l1 penalized negative binomial regressions.  Statistics & Probability Letters 165, 108816.
  * Yang et al. (2021) Yang, X., S. Song, and H. Zhang (2021).  Law of iterated logarithm and model selection consistency for generalized linear models with independent and dependent responses.  Frontiers of Mathematics in China 16(3), 825–856.
  * Yee (2015) Yee, T. W. (2015).  Vector generalized linear and additive models: with an implementation in R.  springer.
  * Zhang and Chen (2021) Zhang, H. and S. X. Chen (2021).  Concentration inequalities for statistical inference.  Communications in Mathematical Research 37(1), 1–85.
  * Zhang and Jia (2022) Zhang, H. and J. Jia (2022).  Elastic-net regularized high-dimensional negative binomial regression: Consistency and weak signals detection.  Statistica Sinica 32, 181–207.
  * Zhang et al. (2018) Zhang, H., K. Tan, and B. Li (2018).  Com-negative binomial distribution: modeling overdispersion and ultrahigh zero-inflated count data.  Frontiers of Mathematics in China 13(4), 967–998.