Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
38,175 characters · 11 sections · 20 citation commands
Heterogeneous Overdispersed Count Data Regressions via Double Penalized Estimations
\def\spacingset#1{ {#1}} \spacingset{1}
\if00 \fi
\if10 {
} \fi
{\it Keywords:} Negative binomial regressions; Heterogeneous count data regression; Estimation of dispersion parameter; Oracle inequalities.
\spacingset{1.8}
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 1996Regression.
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 dai2013maximum,zhang2018negative, 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 xie2020consistency. Recently, there has been much research on the high-dimensional NBR model, such as Qiu17, weissbach2020consistency, tian2020seasonal. 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 Hilbe11. 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 yee2015vector,nguelifack2019robust. 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 wang2016penalized, 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.
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
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:
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
For notation simplicity, we denote
for any measurable function $f$.
Let $\theta=(\theta^{(1)\top},\theta^{(2)\top})^{\top}\in \mathbb{R}^{2p}$,the log-likelihood is
We use the negative log-likelihood as the loss function $\gamma$, and define
Denote $\partial_j := \frac{\partial}{\partial \theta^{(j)}}$, $j = 1, 2$, the score function for $\theta^{(1)}$ is
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
So from now, we will suppose the true parameter is $\theta^*$.
The lasso estimator under our circumstance is defined as
where $\lambda > 0$ is the tuning parameter and the weighted norm is defined by
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 huang2021weighted. Equation ((ref)) 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$.
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^{\dag}=(\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.
The proposed algorithm can do variable selection and dispersion estimation simultaneously. Similar studies on high-dimensional NBR models include wang2016penalized, 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.
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(
\right) \in \mathbb{R}^{2 \times 2p}$, we can assume the value space $\mathcal{S}$ for $s := z \theta$ is bounded and satisfies
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 zhang2022elastic, to solve $\ell_1$-penalization problem. Define the local stochastic Lipschitz constant by
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
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 (adamczak2008tail).
{ {
}}
{ {
It is worthy to note that the $M_{y,n}$ in Corollary (ref) 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 (ref). Corollary (ref) combined Lemma (ref) will give the following corollary as a step more.
}}
Theorem (ref) 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.
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 ((ref)) 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 bickel2009simultaneous, which says that for some $1 \leq s \leq p$ and $K > 0$,
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, bickel2009simultaneous 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 candes2007dantzig as follows
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 ((ref)) can guarantee a good consistent property.
From Theorem (ref), we can get the following corollary regarding the consistency of $\widehat{\theta}$ immediately by the property of the supreme empirical process.
In this section, we evaluate the finite sample performance of the proposed method. The response is generated from the negative binomial regression model ((ref)) with
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 {\bf lbfgs} is required to solve the optimization problem.
{\bf 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 (ref) 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.
{\bf 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 (ref) 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.
In this section, we apply the proposed method to the dataset of German health care demand. The data was employed in riphahn2003incentive and could be downloaded on \url{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 (ref) 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.
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 yang2021law. 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)} \ne 0 . $$ The issues concerning the hypothesis testing are via the debiased Lasso estimator; see shi2019linear 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 xie2021aggregating,cui2021directional.