Law of the Iterated Logarithm and Model Selection Consistency for GLMs with Independent and Dependent Responses

TL;DR

This paper establishes the law of the iterated logarithm for maximum likelihood estimators in generalized linear models with both independent and dependent responses, and applies it to prove the strong consistency of model selection criteria like BIC.

Contribution

It extends the law of the iterated logarithm to GLMs with dependent responses and demonstrates the resulting strong consistency of penalized likelihood model selection methods.

Findings

LIL holds for MLE in GLMs with dependent responses.

Penalized likelihood criteria can consistently select the correct model.

Simulation confirms BIC's selection consistency.

Abstract

We study the law of the iterated logarithm (LIL) for the maximum likelihood estimation of the parameters (as a convex optimization problem) in the generalized linear models with independent or weakly dependent ($\rho$-mixing, $m$-dependent) responses under mild conditions. The LIL is useful to derive the asymptotic bounds for the discrepancy between the empirical process of the log-likelihood function and the true log-likelihood. As the application of the LIL, the strong consistency of some penalized likelihood based model selection criteria can be shown. Under some regularity conditions, the model selection criterion will be helpful to select the simplest correct model almost surely when the penalty term increases with model dimension and the penalty term has an order higher than $O({\rm{loglog}}n)$ but lower than $O(n)$. Simulation studies are implemented to verify the selection consistency of BIC.

Tables

Table 1. Table 1: Simualtion results for variable selection.

Model	Method	Sample Size	Correct	Overfit	Underfit	MSE
NBR ($\beta =0.5,\theta =10$)	BIC	100	0.8920	0.1053	0.0037	0.0164
		300	0.9427	0.0573	0.0000	0.0064
	AIC	100	0.5617	0.4373	0.0010	0.0287
		300	0.5863	0.4137	0.0000	0.0197
Probit model($\beta =0.5$)	BIC	100	0.7106	0.0816	0.2078	0.0190
		300	0.9444	0.0000	0.0556	0.0094
	AIC	100	0.5420	0.4164	0.0416	0.0155
		300	0.5928	0.0000	0.4072	0.0193
Dependent LM (MR(2) error)	BIC	100	0.8460	0.1002	0.0538	0.0178
		300	0.9468	0.0532	0.0000	0.0104
	AIC	100	0.5590	0.4330	0.0080	0.0222
		300	0.5928	0.0000	0.4072	0.0193
Dependent LM (MR(3) error)	BIC	100	0.8366	0.1024	0.0610	0.0181
		300	0.9510	0.0490	0.0000	0.0110
	AIC	100	0.5642	0.4254	0.0181	0.0254
		300	0.5890	0.4110	0.0000	0.0116

Equations

$$\underset{n\to \infty }{\mathrm{limsup}}\frac{\parallel \hat{\beta }-{\beta }_0\parallel }{\sqrt{n^{-1}\log \log n}}=da.s.,$$

$$dF(y_i)=c(y_i)\exp \{y_i{\theta }_i-b({\theta }_i)\}d\mu (y_i),c(y_i)>0,i=1,\dots ,n$$

$${\theta }_i=u(x_i^t{\beta }_0).$$

$$E(y_i|x_i):={\mu }_i=\dot{b}({\theta }_i)=\dot{b}(u(x_i^t{\beta }_0)),i=1,2,\cdots ,n.$$

$$l_n(\beta ):=\sum _{i=1}^n[y_{i}u(x_i^t\beta )-b(u(x_i^t\beta ))].$$

$$S_n(\beta ):=\frac{\partial l_n(\beta )}{\partial \beta }=\sum _{i=1}^n{x}_i\dot{u}(x_i^t\beta )[y_i-\dot{b}(u(x_i^t\beta ))].$$

$$S_n(\beta )=\sum _{i=1}^n{x}_i\dot{u}(x_i^t\beta )[y_i-\dot{b}(u(x_i^t\beta ))]=\sum _{i=1}^n{x}_i[y_i-E(y_i)].$$

$$l(y,\beta )=\frac{1}{n}\sum _{i=1}^n\{y_i\log \Phi (x_i^t\beta )+(1-y_i)\log [1-\Phi (x_i^t\beta )]\}.$$

$$u(x_i^t\beta )=\log [\Phi (x_i^t\beta )/(1-\Phi (x_i^t\beta ))],b(u(x_i^t\beta ))=\log (1-\Phi (x_i^t\beta )),$$

$$f(y_i|\theta ,{\mu }_i)=\frac{\Gamma (\theta +y_i)}{\Gamma (\theta )y_i!}(\frac{{\mu }_i}{\theta +{\mu }_i}{)}^{y_i}(\frac{\theta }{\theta +{\mu }_i}{)}^{\theta },i=1,2,\cdots ,n.$$

$$l(y;\beta )$$

$$u(x_i^t\beta )=x_i^t\beta -\log (\theta +\exp \{x_i^t\beta \}),b(u(x_i^t\beta ))=\theta \log (\theta +\exp \{x_i^t\beta \})$$

$$l_n(\beta ):=\sum _{i=1}^n{w}_i[y_{i}u(x_i^t\beta )-b(u(x_i^t\beta ))],$$

$${\hat{\beta }}_n=\arg \underset{\beta \in {\mathbb{R}}^p}{\max }\{\sum _{i=1}^n{w}_i[y_{i}u(x_i^t\beta )-b(u(x_i^t\beta ))]\}.$$

$$S_n(\beta ):=\frac{\partial l_n(\beta )}{\partial \beta }$$

$$\frac{{\partial }^2{l}_n(\beta )}{\partial \beta \partial {\beta }^t}=\sum _{i=1}^n{w}_i\{x_i\ddot{u}(x_i^t\beta )x_i^t[y_i-\dot{b}(u(x_i^t\beta ))]-x_i{\dot{u}}^2(x_i^t\beta )\ddot{b}(u(x_i^t\beta ))x_i^t\},$$

$$I_n(\beta ):=-E\frac{{\partial }^2{l}_n(\beta )}{\partial \beta \partial {\beta }^t}=-\sum _{i=1}^n{w}_i\{x_i\ddot{u}(x_i^t\beta )x_i^t[E{y}_i-\dot{b}(u(x_i^t\beta ))]-x_i{\dot{u}}^2(x_i^t\beta )\ddot{b}(u(x_i^t\beta ))x_i^t,$$

$$I_n({\beta }_0)=\sum _{i=1}^n{w}_i{\dot{u}}^2(x_i^t{\beta }_0)\ddot{b}(u(x_i^t{\beta }_0))x_i{x}_i^t.$$

$${\mu }_i:=E(y_i|x_i)=\dot{b}(u(x_{i\alpha }^t\beta (\alpha ))),i=1,2,\cdots ,n.$$

$${\Gamma }_c=\{\alpha :{\beta }_{0j}(\alpha )=0,\text{for any}j\notin \alpha \}.$$

$${\Gamma }_w=\{\alpha :{\beta }_{0j}(\alpha )\ne 0,\text{for some}j\notin \alpha \}.$$

$$S_n(\alpha ):$$

$$:$$

$$\hat{\alpha }=\underset{\alpha \in \mathcal{A}}{\mathrm{argmin}}{S}_n(\alpha ).$$

$$P\{\underset{n\to \infty }{\lim }\hat{\alpha }={\alpha }_0\}=1.$$

$$I_n^{w^{\ast }}({\beta }_0)=\sum _{k=1}^n{w}_k^{\ast }{\dot{u}}^2(x_k^t{\beta }_0)\ddot{b}(u(x_k^t{\beta }_0))x_k{x}_k^t.$$

$$c_{l}n{I}_p\prec I_n^{w^{\ast }}({\beta }_0)\prec c_{u}n{I}_p,$$

$$\underset{k\ge 1}{\sup }E|e_k{|}^3=\underset{k\ge 1}{\sup }E|y_k-\dot{b}(u(x_k^t{\beta }_0)){|}^3\propto \underset{k\ge 1}{\sup }[\ddot{b}(u(x_k^t{\beta }_0)){]}^{3/2}<\infty .$$

$$\Vert \hat{\beta }(\alpha )-{\beta }_0(\alpha )\Vert =O(\sqrt{n^{-1}\log \log n})a.s.,$$

$$\underset{n\to \infty }{\mathrm{limsup}}\frac{\Vert \hat{\beta }(\alpha )-{\beta }_0(\alpha )\Vert }{\sqrt{n^{-1}\log \log n}}=ba.s.,$$

$$0\le l_n(\hat{\beta }(\alpha ))-l_n(\beta ({\alpha }_0))=O(\log \log n)a.s.,$$

Full text

Keywords:  generalized linear models; weighted scores method; non-natural link function; model selection; consistency; weakly dependent.

1 INTRODUCTION

Originating from Nelder (1972), generalized linear models (GLMs) are remarkable and synthetic extensions of linear models. GLMs are often classified into two classes in references. The first type is the GLM with a natural link function (canonical link function), such as the binomial regression (logistic regression) and the Poisson regression. GLMs endowed with non-natural link functions become the second type of GLMs, including the probit model and the negative binomial regression, which are more complex to analyze, see Fahrmeir and Kaufmann (1985), Chen (2011). Recently, GLMs have become the very popular regression models in the big data era, see Efron and Hastie (2016).

Under some regularity conditions, asymptotic normality for maximum likelihood estimator (MLE) in GLMs with both natural and non-natural link functions, was established by Fahrmeir and Kaufmann (1985). However, asymptotic normality is a type of convergence which is weaker than the strong limit theorems. A lot of efforts have been done in studying strong consistency in terms of the law of iterated logarithm (LIL). In linear models, Lai and Wei (1982) obtained a general LIL for weighted sums of the independent regression noises with zero means and common variance, which can be applied to least squares estimation. He and Wang (1995) proved the LIL for a general class of M-estimators and then derived the Bahadur representation. Fang (1998) studied the LIL for the MLE of general nonhomogeneous Poisson processes. If there are measurement errors in covariates, under some regularity conditions, Miao and Yang (2011) established the LIL for the least squares estimation in simple linear models.

For many practical regression problems, redundant or irrelevant covariates will come into the models. If those irrelevant covariates are enclosed, the efficiency of estimators will more or less be impaired. It is indispensable to do variable selection after obtaining the regression coefficients in scientific analysis. Usually, when we consider proving the consistency of model selection like the BIC, weak or even strong consistency of the estimated coefficient is not enough if convergence rate cannot be obtained. Thus, it is worth concerning the rate of strong consistency of the estimator as is the LIL. Under i.i.d. noise assumptions (may be non-normal), Rao and Wu (1989) applied the LIL to the selection consistency of linear models with non-normal noise. Similarly, Wu and Zen (1999) studied the Huber’s M-estimator for linear models.

This work focuses on studying a set of the more flexible GLMs where more link functions and exponential family responses are permitted. In some simple GLMs, canonical links do not always provide the best fit. For example, as the non-canonical link GLM, the negative binomial regression can model over-dispersed count data while the Poisson regression (a GLM with the canonical link function) requires equal dispersion, which is inapplicable in practice. Generally, there is no apriori reason to explain why a canonical link should be used, and in many cases a non-canonical link is more suitable (McCullagh and Nelder (1989) and Czado and Munk (2000)). Note that Qian and Wu (2006)(merely binomial regression) only considered some special cases of GLMs. Our results are the extensions of Qian and Wu (2006) but with slightly differences (see Remark1, (H.5) and (H.6)). Meanwhile, we require fewer assumptions on the Fisher information matrix of true regression parameters.

Let $\beta_{0}:=(\beta_{01},\ldots,\beta_{0p})^{t}$ be the true value of regression parameter $\beta$ and $\|\cdot\|$ be the Euclidian norm. In order to get the model selection consistency, we need to show the LIL for the MLE $\hat{\beta}$ of GLMs. Assume that there exists a constant $d>0$ such that

[TABLE]

under some conditions (see Section 3.3 for details). Then we have the convergence rate $\parallel\hat{\beta}-\beta_{0}\parallel=O(\sqrt{n^{-1}\log\log n})~{}a.s..$ Hence the strong consistency of $\hat{\beta}$ is directly implied by the LIL.

Moreover, this work gives a general consideration of model selection consistency in GLMs with non-natural link functions and the responses are allowed to be weakly dependent including $\rho$-mixing and $m$-dependent. Our contributions are to extend previous works under some mild conditions, which lead to desired performance of model selection, while the previous studies about GLMs seldom consider the non-natural link function and weakly dependent responses. Our generalization covers a wide range of GLMs.

This paper is organized as follows. In Section 2, a review of exponential family GLMs is presented and the problem of estimating the coefficients in GLMs by weighted scores equation is discussed. In section 3, under regularity assumptions, we give the LIL for the weighted scores estimates for independent or weakly dependent responses. Checking the model selection consistency is equivalent to evaluating whether the order of penalty term is between $O({\rm{loglog}}n)$ and $O(n)$. For independent or weakly dependent exponential family responses, it is easy to show the model selection consistency by applying the LIL proposed in Section 3. In section 4, we run the simulation for 500 times, and exemplifies that the BIC is consistent while the AIC is not consistent. The detailed proofs of the theorems and lemmas are demonstrated in APPENDIX.

2 GENERALIZED LINEAR MODELS

We expect readers to be familiar with exponential families as prerequisites. Here we give a brief review that makes our exposition self-contained.

Let $\{(y_{i},x_{i})\}_{i=1}^{n}$ be $n$ independent sample data pairs. Let ${x_{i}}^{,}s\in\mathbb{R}^{p}$ be fixed covariates and assume that the responses ${y_{i}}^{,}s$ follow the density of the exponential families

[TABLE]

where ${\theta_{i}}\in\Theta=\{\eta:\int c(y)\exp\{y\eta\}\mu(dy)<\infty\}$, $c(y_{i})$ is free of $\theta_{i}$ and $\mu(\cdot)$ is the dominated measure.

Let $\dot{b}(\cdot)$ and $\ddot{b}(\cdot)$ be the first and second derivatives of $b(\cdot)$, respectively. A notable property of the expectation and variance of ${y_{i}}^{,}s$ in (1) is $E(y_{i})=\dot{b}(\theta_{i})$ and $Var(y_{i})=\ddot{b}(\theta_{i})$. Note that the first and second derivatives of $b(\cdot)$ are bounded in any compact set $K\subset\Theta$. For more details of exponential families, please see Brown (1986).

For GLMs, we assume that the effect of covariates $x_{i}$ on responses $y_{i}$ can be observed by the corresponding regression coefficient $\beta$, and the relationship between $\theta_{i}$ and $x_{i}^{t}\beta$ is expressed by the following relation function

[TABLE]

Therefore, $\theta_{i}$ is determined by the unknown $\beta$ and the given $x_{i}$. Here we suppose that $x_{i}$’s are non-random vectors of dimension $p$. The $\beta\in\mathbb{R}^{p}$ is called regression coefficient which will be estimated.

Further, we have

[TABLE]

Let $g$ be the link function such that $g(\mu_{i})=x_{i}^{t}\beta_{0}$. Since $\mu_{i}=g^{-1}(x_{i}^{t}\beta_{0})$ with $\mu_{i}=\dot{b}(u(x_{i}^{t}\beta_{0}))$, we immediately acquire the expression $u(t)=\dot{b}^{-1}(g^{-1}(t))$.

The likelihood function of $\{(y_{i},x_{i})\}_{i=1}^{n}$ is the product of $n$ terms in (1). Taking the logarithm, we get the log-likelihood function

[TABLE]

Then the score function is defined by

[TABLE]

For more discussion of GLMs, readers are suggested to see McCullagh and Nelder (1989), Fahrmeir and Tutz (2001), Chen (2011).

2.1 Examples of Non-natural Link

For GLMs with the natural link function $g(t)=\dot{b}^{-1}(t)$, we have the identity function $u(t)=t$. Thus the first term in (2) is a linear function of $x_{i}^{t}\beta$ and the score function is

[TABLE]

Furthermore, we illustrate some GLMs with non-natural link functions.

1. Probit Model

Let $\Phi(t)$ be the cumulative distribution function of $N(0,1)$. Suppose the $n$ independent observed data sets $\{y_{i},x_{i}\}_{i=1}^{n}$ are from a logit model, and $y_{i}$’s are independent Bernoulli distributed with $P(y_{i}=1|x_{i};\beta)={\Phi(x_{i}^{t}\beta)},~{}i=1,2,...,n$. The average log-likelihood function for data $\{y_{i},x_{i}\}_{i=1}^{n}$ is

[TABLE]

Like the logistic regression, the $u(\cdot)$ and $b(\cdot)$ in the probit model follow parametric equations:

[TABLE]

where $\Phi(t)$ is the cumulative distribution function of the standard normal distribution.

2. Negative Binomial Regression

It is common that negative-binomial GLMs are more plausible than Poisson regression for modelling overdispersion count data. Negative binomial regression (NBR) assumes that the overdispersed response data are modelled by two-parameter distribution:

[TABLE]

where $\theta$ is known (can be estimated previously), $E(y_{i}|x_{i})=\mu_{i}$, $Var(y_{i}|x_{i})=\mu_{i}+\frac{\mu_{i}^{2}}{\theta}$.

Suppose that the relationship between the mean parameter and covariates is given by $\log(\mu_{i})=x_{i}^{t}\beta$. Then the logarithm of the maximum likelihood function for NBR is

[TABLE]

Thus the connection of $u(\cdot)$ and $b(\cdot)$ for NBR is

[TABLE]

where $\mu_{i}=\exp\{x_{i}^{t}\beta\}$.

NBR plays an important role in modern applications relating to count data, see Zhang and Jia (2017) and references therein.

2.2 Weighted scores equations

When estimating regression coefficients, a well-known robust approach is the weighted-likelihood method which assigns a sequence of weights to perturb the contribution of each sample in the log-likelihood function. Based on some weight functions (see Markatou et al. (1998)), the simultaneous changes in the weights of $n$ samples producing the robust estimators for GLMs. It is a method of mitigating the influence of outliers (or say the leverage points, see Pages 59 of van der Vaart (1998)).

If the log-likelihood (2) is replaced by the following weighted log-likelihood

[TABLE]

the maximum likelihood estimation for the weighted-likelihood is

[TABLE]

Here $\{w_{i}\}_{i=1}^{n}$ are any uniformly bounded positive weights that are independent of response $y_{i}$ for all $i$ (The weights may depend on $\{x_{i}\}_{i=1}^{n}$).

Then the weighted MLE is a robust estimation for some suitable weights $\{w_{i}\}_{i=1}^{n}$. If $w_{i}\equiv$ constant for all $i$, then the unweighted MLE is exactly the common GLMs, which may not have the robust property.

Consider vector derivative with respect to $\beta$ by letting $\frac{\partial l_{n}(\beta)}{\partial\beta}=\{\frac{\partial l_{n}(\beta)}{\partial\beta_{1}},\cdots,\frac{\partial l_{n}(\beta)}{\partial\beta_{p}}\}^{t}$. Then weighted score is given by

[TABLE]

In particular, we have $S_{n}(\beta_{0})=\sum\limits_{i=1}^{n}w_{i}x_{i}\dot{u}(x_{i}^{t}\beta_{0})[y_{i}-E(y_{i})]$.

From $\frac{\partial l_{n}(\beta)}{\partial\beta}$, assume that $\hat{\beta}$ is a unique solution by the weighted scores equation $S_{n}(\beta)=0$. Note that the Hessian matrix of $\beta$ is

[TABLE]

and the Fisher information at $\beta$ is

[TABLE]

where the last equality is due to $E(y_{i})=\dot{b}(u(x_{i}^{t}\beta))$. In particular, we have

[TABLE]

It is easy to see that the Fisher information $I_{n}(\beta)$ is semi-positive definite. Then, $\hat{\beta}$ makes the likelihood function get the maximum.

3 MODEL SELECTION CONSISTENCY

3.1 Preliminaries and Notations

In order to examine the predictability of the model, in this section we discuss the model selection methods by including or excluding variables. To be specific, for $i=1,\cdots,n$, the question is “Which optimal subset of $x=({x_{i1}},\cdots,{x_{ip}})^{t}$ will enter into the regression function $E(y_{i}|x_{i})$ by some model selection approaches?”. In the existing literatures, there are several commonly used model selection criteria such as AIC (Akaike Information Criterion, Akaike (1973)), BIC (Bayesian Information Criterion, Schwarz (1978)) and SCC (Stochastic Complexity Criterion, Rissanen (1989)), whose theoretical background is rooted in the information theory or the Bayesian analysis.

To answer the question, we now build a mathematical framework to get optimal sub-model. Let $\alpha$ be a subset of $\{1,2,\ldots,p\}$ and define $\beta(\alpha)$ (or $x_{\alpha}$) to be the sub-vector of $\beta$ (or $x$) indexed by the integers in $\alpha$. The dimension of $\beta({\alpha})$ is denoted by $p_{\alpha}$. Therefore, for simplicity, we say that the sub-model corresponding to $\alpha$ is called the $\alpha$ sub-model or the candidate model. It is given by the mean of the predictor

[TABLE]

Note that $\alpha$ sub-model is not necessarily a correct model in a manner that $E({y_{i}}|{x_{i}})$ is not constantly and precisely equal to $\dot{b}(u(x_{i\alpha}^{t}\beta(\alpha)))$. Let ${\cal A}$ be the collection of all subsets of ${\{1,2,\cdots,p\}}$. Then there are $|{\cal A}|=2^{p}$ candidate models for screening. If $\beta_{0i}=0$ for all i not in alpha then model alpha is called "correct". The set of all correct models is denoted $\Gamma_{c}$. The set $\alpha_{0}=\{i:\beta_{0i}\neq 0\}$ is the smallest correct model. The set of all models which are not correct is denoted $\Gamma_{w}$. Then model $\alpha$ is called a correct model, denoted by

[TABLE]

$\Gamma_{c}$ is the set of all correct models. However, there could be more than one correct models. Many models belonging to $\Gamma_{c}$ may also contain some redundant variables $x_{ij}$ that have no effect on response $y_{i}$. The remaining candidate models can be collected into

[TABLE]

Here $\Gamma_{w}$ is the set of all wrong models, each of which misses at least one non-zero variable $x_{ij}$ that has an effect on $y_{i}$.

Therefore, it is quite simple to obtain an objective function based on the GLMs likelihood, and use any of the above model selection methods to evaluate the accuracy of all sub-models, unless $p$ can be large enough that implementation of the model selection is computationally prohibitive. In this work, we only study the fixed dimension case, while consistent model selection criteria on the increasing- or high-dimensional case have been studied by Kim and Jeon (2016). Now define

[TABLE]

where each $\alpha$ sub-model is estimated via the MLE from the analogue of (4) for the subvector $\beta(\alpha)$.

Here, the first term in (3.1) is the negative log-likelihood that is used to detect the superiority of the model. The second term is a penalty term used to measure the complexity of the model. We assume that $C(n,\hat{\beta}(\alpha))$ is an increasing function of the sub-model dimension ${p_{\alpha}}$. For AIC, $C(n,{{\hat{\beta}}_{n}}(\alpha))={p_{\alpha}}$. For BIC, $C(n,{{\hat{\beta}}_{n}}(\alpha))=\frac{1}{2}{p_{\alpha}}\log n$. For SCC, $C(n,\hat{\beta}_{n}(\alpha))=\log|I_{n}(\hat{\beta}_{n}(\alpha))|/2+\sum\limits_{i=2}^{p_{\alpha}}\log(|\hat{\beta}_{n}(\alpha)_{i}|+\varepsilon n^{-1/4})$ ,where $I_{n}(\beta(\alpha))$ is the Fisher information matrix of $\beta(\alpha)$ and $\varepsilon$ is a specific value to ensure the invariance of the model selection criterion.

Minimize (3.1) of all candidate models to get the optimal model, i.e.,

[TABLE]

Under the above criteria, the better the model is fitted, the less complex the model is, i.e., the smaller ${S_{n}}(\alpha)$ is, the better the model selection effect is. The penalized criterion in (7) shows that the model with good predictability should not only enjoy a small fitting deviation (by the first term in (3.1)), but also not incorporate unessential variables, i.e., $C(n,\hat{\beta}_{n}(n))$ should be quite small. We know that a model with the smallest ${S_{n}}(\alpha)$ among all candidate models would be optimal.

Let $\alpha_{0}$ be the correct model with the smallest size, which is called the simplest correct model. That is, the simplest correct model includes exactly all nonzero components of $\beta({\alpha_{0}})$. For the simple presentation, we assume that the simplest correct model is unique. Let ${\hat{\alpha}}$ be the estimate from the model selection criterion (7), then the model selection procedure is said to be strongly consistent if it can select the simplest correct model almost surely, say

[TABLE]

The main goal of Section 3 is to evaluate the consistency of the model selection methods that can select the optimal model from all candidate models. Later in the article, some regularity conditions are used to establish model selection criteria and to give some asymptotic results.

3.2 Regularity Conditions

We give some general assumptions in advance.

• •

(H.1): Let ${\lambda_{1}}\{S\}\leq\cdots\leq{\lambda_{p}}\{S\}$ be the ordered $p$ eigenvalues of a $p\times p$ symmetric matrix $S$. Then

(i) Let $\{I_{n}(\beta_{0})\}$ be Fisher information of $\beta$ given by (2.2). We have $\lim\limits_{n\rightarrow\infty}\lambda_{j}\{I_{n}(\beta_{0})\}=\infty$, $j=1,\ldots,p$; (ii) There exist positive constants $d_{1}$,$d_{2}$, that satisfy $d_{1}n\leq\lambda_{p}\{I_{n}(\beta_{0})\}\leq d_{2}n$.

• •

(H.2): Bounded elements of the non-random design: For all $i$, we assume that ${\left\|x_{i}\right\|_{\infty}}:=\mathop{\max}\limits_{i,j}\left|{{x_{ij}}}\right|\leq L$, where $L$ is a positive constant.

• •

(H.3): We assume that all the coefficients $\beta$ are in parameter space belonging to the space in $\mathbb{R}^{p}$ such that $\mathop{\sup}\limits_{k}\left|{\ddot{u}(x_{k}^{t}{\beta})}\right|,\mathop{\sup}\limits_{k}\left|{\dot{u}(x_{k}^{t}{\beta})}\right|,\mathop{\sup}\limits_{k}\ddot{b}(u(x_{k}^{t}{\beta}))<\infty$.

• •

(H.4): In the weighted GLMs, if the weights are replaced by some uniformly bounded new weights ${w_{k}^{*}}$, $k=1,\cdots,n,$ which are all uniformly bounded, then there exists a positive constant $W$ that satisfies $\mathop{\max}\limits_{1\leq k\leq n}w_{k}^{*}\leq W$ for all $n$. For notation simplicity, we denote the Fisher information of $\beta$ with old weights in (2.2) as $I_{n}^{w}(\beta):={I_{n}}(\beta)$. We define the weighted Gram matrix as

[TABLE]

The notation $A\prec B$ means that $B-A$ is non-negative definite. The following restricted eigen-value condition is true:

[TABLE]

where $I_{p}$ is the $p$ dimensional identity matrix, and $c_{l}$, $c_{u}$ are positive constants.

Conditions (H.1) below are also proposed in Wu and Zen (1999), Qian and Wu (2006) for some special cases of GLMs. Moreover, the combination of (H.1) and (H.2) renders a reasonable convergence rate of estimated coefficients. The bounded regressors assumption (H.2) is common in some GLM references, see page 46 of Fahrmeir and Tutz (2001) and Section 2.2.7 of Chen (2011) as examples. When we preprocess the raw covariates, “zero mean” and “one variance” standardizations are required, which evidently and approximately imply the boundedness assumption of covariates. If some predictors have heavy tailed distributions which may be collected in economics and finance, we could take the some-transform such as $f(x)=c\frac{\exp(x)}{1+\exp(x)}$ (where $c$ is a constant) and thus the transformed predictors are approximately seen as bounded variables. For the random design case, we can get all the results by conditioning on $X$. The (H.3) is important for the consistency results (see Shao (2003)) since the gradient and Hessian for the weighted likelihood in (4) should be bounded to make sure that the optimization problem has a good solution. In the proof, we need to ensure that all the $\mathop{\sup}\limits_{k}\left|{\ddot{u}(x_{k}^{t}{\beta_{0}})}\right|,\mathop{\sup}\limits_{k}\left|{\ddot{u}(x_{k}^{t}{\beta_{1}})}\right|,\mathop{\sup}\limits_{k}|{\dot{u}(x_{k}^{t}{\beta_{0}})}|,\mathop{\sup}\limits_{k}|{\dot{u}(x_{k}^{t}{\tilde{\beta}_{1}})}|,\mathop{\sup}\limits_{k}|{\dot{u}(x_{k}^{t}{\breve{\beta}_{1}})}|$, and $\mathop{\sup}\limits_{k}\ddot{b}(u(x_{k}^{t}{\beta_{0}}))$ are finite, where the values of $\beta_{1}$ and $\tilde{\beta}_{1}$ are between the line of $\beta$ and $\beta_{0}$, and $\beta_{1}$ is not necessarily equal to $\tilde{\beta}_{1}$. In addition, the value of $\breve{\beta}_{1}$ is between the line of $\hat{\beta}$ and $\beta_{0}$. As for (H.4), with (H.3) it is the theoretical guarantee relating to the proofs of the LIL of the estimator.

Remark 1**.**

By (H.3) there is a compact subset $K$, not depending on $n$, of the interior of $\Theta$ such that $u(x_{k}^{t}{\beta_{0}})$ belongs to $K$ for all $k$. Let ${e_{k}}={y_{k}}-\dot{b}(u(x_{k}^{t}{\beta_{0}}))$, the Lemma 6.1 in Rigollet (2012) implies

[TABLE]

The moments condition (8) is crucial to explain the finite moments behavior of the exponential families. This condition means, in most cases, that there is a compact subset $K$, not depending on $n$, of the interior of $\Theta$ such that $u(x_{k}^{t}{\beta_{0}})$ belongs to $K$ for all $k$. Let $\eta>0$. Condition $\mathop{\sup}\limits_{k\geq 1}E|{e_{k}}{|^{2+\eta}}<\infty$ is used in Yin et al. (2006) for establishing the strongly consistent maximum quasi-likelihood estimates.

3.3 Independent Observations

First, we give the following useful strong limit theorems.

Theorem 3.1**.**

(Law of the Iterated Logarithm) Assuming conditions (H.1) to (H.4) are satisfied, for any correct model $\alpha\in\Gamma_{c}$,

[TABLE]

Furthermore, there is a positive constant $b>0$ such that for arbitrary $\alpha\in\Gamma_{c}$

[TABLE]

that is, the MLE $\hat{\beta}(\alpha)$ satisfies the LIL.

Once Theorem 3.1 has been constructed, the following almost surely discrepancy between estimated correct sub-model (or incorrect sub-model) log-likelihood and true log-likelihood can be derived as well.

Theorem 3.2**.**

Given the same conditions as Theorem 3.1, for any correct model $\alpha\in\Gamma_{c}$,

[TABLE]

where ${l_{n}}(\hat{\beta}(\alpha))=-\sum\limits_{i=1}^{n}{{w_{i}}}[b(u(x_{i\alpha}^{t}\hat{\beta}(\alpha)))-{y_{i}}u(x_{i\alpha}^{t}\hat{\beta}(\alpha))]$.

From the conclusion of Theorem 3.2, we can see that the maximum log-likelihood for any correct model $\alpha\in\Gamma_{c}$ is almost surely greater than the unknown log-likelihood of the true model, and the bounds of their differences are almost surely limited by $|O(\log\log n)|$.

Theorem 3.3**.**

Under the same conditions as Theorem 3.1, for any incorrect model $\alpha\in\Gamma_{\omega}$, we have

[TABLE]

where ${l_{n}}({\beta_{0}})=-\sum\limits_{i=1}^{n}{{w_{i}}}[b(u(x_{i}^{t}{\beta_{0}}))-{y_{i}}u(x_{i}^{t}{\beta_{0}})]$.

From Theorem 3.3 we can see that the maximum log-likelihood for all non-correct models $\alpha\in\Gamma_{\omega}$ is almost certainly less than the unknown log-likelihood of the true model. At the same time, when $n$ is sufficiently large, the bounds of their differences will be almost surely limited by $|O(n)|$.

Theorem 3.4**.**

For GLMs with weigted score (2.2), under the conditions of Theorem 3.1, both BIC and SCC selection criteria are strongly consistent, but AIC selection criterion is not strongly consistent.

From Theorem 3.4 derived by Theorem 3.2 and Theorem 3.3, we know that if we use the penalty log-likelihood method to select the model based on the criterion (7) and the penalty term $C(n,\hat{\beta}(\alpha))$ has an order between $|O(\log\log n)|$ and $|O(n)|$. Moreover, the penalty term $C(n,\hat{\beta}(\alpha))$ is the increasing function of the model dimension, so it will almost surely select the simplest correct model that belongs to $\Gamma_{c}$. The phenomenon that the model selection criteria almost surely choose the simplest correct model is called strong consistency of model selection.

3.4 Weakly Dependent Observations

In economics and finance, when responses are collected as time series data with short-term auto-correlation (section 15 in Hansen (2018)), heavily relying on the assumption of independent errors are not logical or reasonable. Whereas, most existing references of studying GLMs (includes linear model) are addressed only for independent data, and dependent data are largely not covered. Some work has paid attention to linear models with correlated errors, see Fan et al. (2016). However, in contrast to linear models with dependent responses, the GLMs emphasizing on dependent errors or dependent responses have only been scarcely investigated, see Kroll (2019). This section focuses on the GLMs endowed with weakly dependent responses, i.e.,

[TABLE]

where ${\varepsilon_{i}}$ are weakly dependent error sequence with zero mean and $\{{x_{i}}\}$ are fixed.

Although, the weakly dependence cannot be solved by MLE approach due to the dependent structure, the likelihood does not enjoy the form of production. However, by applying quasi-likelihood method, we can still be able to get a desired estimator which is consistent under some regularity conditions. The quasi-likelihood estimator $\hat{\beta}_{\rm Q}(\alpha)$ is the solution to follow estimating equation:

[TABLE]

The inference from the solution of (13) is a typical example and is the maximum likelihood estimate for independent responses. More details about quasi-likelihood can be seen in Fahrmeir and Tutz (2001), Chen (2011).

Next, we still want to apply LIL that would imply model selection consistency by using the similar approach mentioned in the case for independent observations. In order to present the dependence structure, we point out some notations and definitions for measuring dependence.

Let $(\Omega,\mathcal{F},P)$ be a probability space, and $\mathcal{B}$, $\mathcal{C}$ two sub-$\sigma$-fields of $\mathcal{F}$. Let ${L_{2}}({\cal B})$ be a set of all ${\cal B}$-measurable random variables with finite 2nd moments.

We introduce the notation of $\rho$-mixing and $\alpha$-mixing coefficients (strong mixing coefficient):

[TABLE]

Let $\mathcal{F}_{a}^{b}=\sigma({X_{t}},a\leq t\leq b)$ and $\mathbb{N}$ be the set of all positive integers. The random sequence $\{{X_{t}},t\in\mathbb{N}\}$ is said to be $\alpha$-mixing or strongly mixing, if

[TABLE]

Next, a sequence of random variables $\{{X_{t}},t\in\mathbb{N}\}$ is called a $\rho$-mixing process if

[TABLE]

A sequence of strictly stationary random variables $\{{X_{t}},t\in\mathbb{N}\}$ are said to be $m$-dependent, if for any two subsets $B,C\subset\mathbb{N}$ and $\mathop{\inf}\limits_{{t_{1}}\in B,{t_{2}}\in C}\left|{{t_{1}}-{t_{2}}}\right|>m$, the process $\{{X_{t_{1}}},{t_{1}}\in B\}$ and $\{{X_{t_{2}}},{t_{2}}\in C\}$ are independent.

For detailed theories and examples (such linear time series) about $\alpha$-mixing process, $\rho$-mixing process and other mixing processes, we refer reading to Lin and Lu (1997) and Bosq (1998). Notice that, we have $\alpha({\cal A},{\cal B})\leq\frac{{\rm{1}}}{{\rm{4}}}\rho({\cal A},{\cal B})$ and then the $\rho$-mixing process implies the $\alpha$-mixing process. For simplicity, we would restrict our study to strictly stationary sequences in this section.

To prepare for the asymptotic theories for the weakly dependent case, we need some additional regularity assumptions for our technique proofs.

• •

(H.5): The ${y_{i}}$ satisfies the $\rho$-mixing condition with geometric decay: $\rho(n)=O({r^{-n}})$.

• •

(H.6): The ${y_{i}}$ is $m$-dependent.

Since weak dependency makes the problem more complex, we assume that the design matrix $X$ is viewed as a non-random matrix. A similar assumption is given by Fan et al. (2016). Differ from Condition 3 in Fan et al. (2016), they assumed that the $\alpha$-mixing condition is by exponential decay $\alpha(m)=O({e^{-a{m^{b}}}})$ where $a,b$ are positive constant. Our geometric rate of decay is not sharper than exponential decay, which means that the dependence of our response is allowed to be stronger. Property (H.5) or (H.6) holds for the $\varepsilon_{i}$ if and only if it holds for the $y_{i}$. The same assertion holds for (H.6). It should be note that for (H.6) if an $m$-dependent sequence has $\rho(n)=O({r^{-n}})$ for $n>m$, then the sequence is $\rho$ mixing. But it may not be $\rho$-mixing with geometric decay $\rho(n)=O({r^{-n}})$. The (H.6) implies the $\rho$-mixing with truncated decay: $\rho(n)=O({a_{n}}\cdot 1\{n\leq m\})$ for some sequence $\{{a_{n}}\}$.

Now, we present the main results below which are the same as situation of independent observations.

Theorem 3.5**.**

(Law of the Iterated Logarithm) Assuming conditions (H.1) to (H.4) are satisfied. For weakly dependent GLMs with estimating equation (13), the $\rho$-mixing responses $\{{y_{i}}\}_{i=1}^{n}$ with additional requirement (H.5) or the $m$-dependent responses $\{{y_{i}}\}_{i=1}^{n}$ with additional requirement (H.6), we have

[TABLE]

for any correct model $\alpha\in\Gamma_{c}$.

Furthermore, there is a positive constant $b>0$ such that for any arbitrary $\alpha\in\Gamma_{c}$

[TABLE]

where MLE estimates $\hat{\beta}(\alpha)$ satisfy the law of iterated logarithm.

Once the LIL for weakly dependent GLMs is verified, the following strong convergence for the BIC or SCC model selection is similarly to be obtained, by checking the Theorem 3.2 and Theorem 3.3.

Theorem 3.6**.**

In the weakly dependent GLMs with estimating equation (13), considering the $\alpha$-mixing responses or $m$-dependent $\{{y_{i}}\}_{i=1}^{n}$ error sequences. Under the same conditions of Theorem 3.5, both BIC and SCC criteria are strongly convergent for the model selection, but not strongly convergent for the AIC criterion.

4 SIMULATION STUDY

The purpose for this section is to examine the difference between the performance of BIC and AIC in the variable selection. The BIC criterion for MLE $\hat{\beta}^{nb}(\alpha)$ of NBR is defined as

[TABLE]

where $\hat{\rm{df}}(\hat{\beta}^{nb}(\alpha)):=||\hat{\beta}^{nb}(\alpha)|{|_{0}}$ is the number of coefficients in the model. Similarly, the BIC criterion for MLE $\hat{\beta}^{pb}(\alpha)$ of probit model is given by

[TABLE]

where $\Phi(t)$ is the cumulative distribution function of $N(0,1)$.

For weakly dependent responses, we consider dependent linear model with $MR(p)$ time series ($p=2,3$, see section 15.4 Hansen (2018)) as the error sequence $\{{\varepsilon_{i}}\}$, that is

[TABLE]

We solve (16) (denote $\hat{\beta}^{dl}(\alpha)$) by the OLS method which is equivalent to the estimation from quasi-likelihood. The BIC criterion for dependent linear model is

[TABLE]

At the same time, for comparison, we consider the AIC criterion which replaces the BIC penalty function $\frac{{\log n}}{n}\hat{\rm{df}}[\hat{\beta}(\alpha)]$ by the penalty function $\frac{1}{n}\hat{\rm{df}}[\hat{\beta}(\alpha)]$.

In the simulation, we add three more redundant variables $X_{4i},X_{5i},X_{6i}$ in each observation, which are also independent and distributed $U(0,1)$. The negative binomial regression model and the probit regression model are considered in our simulation and we set the $\mathbf{\beta}$ in the negative binomial regression and the probit model to be 0.5. The $\theta$ in the NBR is 10 in the simulation. As for dependent linear model, we consider the error term which satisfies $MR(2)$ and $MR(3)$ with parameters (0.5, 0.3) and (0.5, 0.3, 0.2).

We run the simulation for 500 times. Each time all ($2^{6}-1=63$) kinds of combinations of the variables denote a model selection to be “correct” if the method successfully selects exactly $\{X_{i1},X_{i2},X_{i3}\}$. If the criterion picks other variables besides $\{X_{i1},X_{i2},X_{i3}\}$, we say it is an “overfit”. The rest is denoted as “underfit”. The MSE is also calculated in our simulation. From table 1 we can see that, in comparison with AIC, BIC criterion always performs better in selecting the true model. At the same time, with the growth of sample size $N$, the “correct” of BIC criterion gradually converges to 1.

Since the AIC/BIC analysis of GLMs have been applied widely to real data, we ignore the real data analysis in this work, see Tutz (2011).

5 DISCUSSION

This paper presents an unified theory for studying strong limit theorems for independent and dependent GLMs responses and we obtain the LIL for the maximum likelihood estimate of the regression coefficients. The convergence rate of regression coefficients estimators is shown to be $O(\sqrt{n^{-1}\log\log n})$ almost surely. This LIL is employed to establish the strong consistency of BIC and SCC selection criteria.

Using the same techniques in this paper, we can also prove that the BIC procedure for sub-sampling based GLMs estimates in Ai et al. (2020) also enjoys the strong consistency of model selection. Further research of this work may be the LIL for the maximum likelihood estimator in errors-in-variables GLMs under mild conditions, where Miao and Yang (2011) only considered the case of simple linear regression.

6 Acknowledgements

The authors would to show their hearty thanks to Professor YueHua Wu for her helpful comments and mentioned references Wu and Zen (1999) and Rao and Wu (1989) which substantially improve this work.

7 APPENDIX

The key of our proofs is the mathematical analysis of the negative likelihood function, and the local quadratic approximation is frequently used. The main difficulty lies in how to standardize the score function and how to establish a consistent bound on the the difference of two log-likelihood functions almost exactly. The techniques of using some concentration inequalities and the error bound of log-likelihood ratio are broadly employed in literatures, see Rao and Wu (1989), Rao and Zhao (1992), Wu and Zen (1999), compare Qian and Wu (2006) for similar considerations.

First, we define a positive sequence $\{\tau_{n}\}$ such that

[TABLE]

In the sequel, we introduce a sequence of $L_{2}$ ball with radius being proportional to $\{\tau_{n}\}$:

[TABLE]

and denote $\partial B_{n}=\{\beta:\|\beta-\beta_{0}\|=\tau_{n}(n^{-1}\log\log n)^{\frac{1}{2}}\}$ as its boundary. Obviously, we have $B_{1}\supset B_{2}\supset B_{3}\supset\cdots\supset B_{n}$.

Further, we define the log-likelihood ratio of two log-likelihoods in terms of $\beta$ and $\beta_{0}$:

[TABLE]

Let $K(t,s)=b(t)-b(s)-\dot{b}(s)(t-s)$, then

[TABLE]

So ${K_{n}}(\beta,{\beta_{0}})$ can be rewritten as

[TABLE]

Our main idea is to assume that the $u(\cdot)$ function is sufficiently smooth, then we can use Taylor expansion for $u(t)$ at $t=s$:

[TABLE]

where $\zeta$ is between $t$ and $s$.

Let $t=x_{k}^{t}\beta$, $s=x_{k}^{t}\beta_{0}$, $\zeta=x_{k}^{t}\beta_{1}$. The value of $\beta_{1}$ is between $\beta$ and $\beta_{0}$. Then by Taylor’s expansion we have

[TABLE]

where

[TABLE]

Appendix A Lemmas

Before presenting the proof of the main results, we give some lemmas which will be useful in the following sections.

The following lemmas that we will apply is a result in convex optimizations.

Lemma A.1**.**

(Lemma 3.3.10 of Wright (2017)) Assume that function $f(\cdot)$ is strongly convex and $\dot{f}(\cdot)$ is Lipschitz continuously differentiable, then $f(x,y)$ has the following upper and lower bounds for arbitrary vectors x, y,

[TABLE]

where $C_{u}$, $C_{l}$ are positive constants.

The following Petrov’s law of iterated logarithm for sums of independent random variables can be found in p274 of Stout (1974).

Lemma A.2**.**

*(Petrov’s LIL) Let $\left(X_{i}\right)_{i\geqslant 1}$ be a sequence of centered independent random variables and define ${S_{n}}=\sum\nolimits_{i=1}^{n}{{X_{i}}}$. Suppose that:

*(i) $\mathop{\lim}\limits_{n\to\infty}\frac{1}{n}{\rm{Var}}{S_{n}}={\sigma^{2}}<\infty$; (ii) $\mathop{\sup}\limits_{i\geq 1}E|X_{i}|^{2+\eta}<\infty$ for some positive $\eta$.

Then*

[TABLE]

Corollary A.1**.**

Under (H.1)-(H.4), we have for $j=1,2,\cdots,p.$

[TABLE]

where $w_{k}$ is the $k$th component of weight $w$, $x_{kj}$ is the $j$th element of vector $x_{k}$, and $I_{n}(\beta_{0})(j,j)$ is the $(j,j)$th element of the Fisher information matrix $I_{n}(\beta_{0})$. The above formula shows that $\{w_{k}\dot{u}(x_{k}^{t}\beta_{0})(y_{k}-\dot{b}(u(x_{k}^{t}\beta_{0})))\}$ meets the LIL, so we have

[TABLE]

where ${1_{p}}:={(1,\cdots,1)^{t}}$.

Corollary A.2**.**

Under (H.1)-(H.4), we have

[TABLE]

Appendix B Proof of Corollary A.1

In the following, we only need to verify (20). Note that the mean and variance of $\{y_{k}\}_{k=1}^{n}$ are $E(y_{k})=\dot{b}(u(x_{k}^{t}\beta_{0})),Var(y_{k})=\ddot{b}(u(x_{k}^{t}\beta_{0}))$ and the definition of $I_{n}(\beta_{0})$, for $j=1,\cdots,p$, let

[TABLE]

According to condition (H.4) and (H.1)(i), we have

[TABLE]

Thus condition (i) of Lemma A.2 is verified.

To check (ii) of Lemma A.2, we only need to use Remark 1 in the Section 3.2.

Appendix C Proof of Corollary A.2

The result (22) comes from (20) and condition (H.1). The proof is similar to Corollary A.1.

Appendix D Proof of Theorem 3.1

Judging from the lemmas above, now we are able to prove Theorem 3.1. For notation simplicity, we write $\beta(\alpha),{x_{k\alpha}}$ as $\beta,{x_{k}}$, respectively.

Firstly, we give the proof of the non-canonical link GLMs.

Let $t=u(x_{k}^{t}\beta)$, $s=u(x_{k}^{t}\beta_{0})$, ${\tilde{\beta}_{1}}$ be the value between $\beta_{0}$ and $\beta$. Then it can be seen from Lemma A.1 that

[TABLE]

where the last equality is from the Taylor’s expansion.

Then

[TABLE]

Denote a weighted Gram matrix as

[TABLE]

Using (H.3) and (H.4), we have ${c_{l}}nI_{p}\prec I_{n}^{{w^{*}}}({\beta_{0}})$ for large $n$. Therefore,

[TABLE]

By $\partial B_{n}=\{\beta:\|\beta-\beta_{0}\|=\tau_{n}\sqrt{n^{-1}\log\log n}\}$, then

[TABLE]

For ${R_{2n}}(\beta)$, according to the Cauchy’s inequality and the LIL, we have

[TABLE]

Then as for ${R_{3n}}(\beta)$, by the Cauchy’s inequality again, it gives

[TABLE]

By Corollary A.2, we get

[TABLE]

So the last equality in (30) is obtained.

Therefore,

[TABLE]

So, for large $n$, we get

[TABLE]

since $\tau_{n}^{2}\log\log n>>\tau_{n}\log\log n$ and $\tau_{n}\to\infty~{}(\text{or~{}some~{}big~{}constant})$.

From Lemma A.1, we see that ${K_{n}}(\beta,{\beta_{0}})$ is convex and then ${K_{n}}(\beta_{0},{\beta_{0}})=0$. The (D) shows that the ${\hat{\beta}}$ minimizes $\frac{1}{n}{K_{n}}(\beta,{\beta_{0}})$ to be contained in a subset of $B_{n}$ almost surely. Then

[TABLE]

Since the sequence $\{\tau_{n}\}$ is chosen to be as slow as possible and tends to infinity (or tends to a constant $C$ such that $O(1)\tau_{n}^{2}\log\log n+O(1)\tau_{n}\log\log n+O(1)n^{-0.5}\tau_{n}^{2}{(\log\log n)^{1.5}}\geq{C^{2}}{\log\log n}$ and $C$ is depended by the constants on (28), (29) and (32), then (34) can be used to get the establishment of (9).

Furthermore, we need to prove

[TABLE]

The following is obtained by contradiction. If it is assumed that the above formula is not established, then

[TABLE]

By applying the result of Lemma A.1, let $\breve{\beta}_{1}$ be the value between $\beta_{0}$ and $\hat{\beta}$, we have

[TABLE]

where the last equality stems from (36).

Similarly, the order of ${R_{2n}}(\hat{\beta}),{R_{3n}}(\hat{\beta})$ can also be drawn

[TABLE]

So under the assumption of (36),

[TABLE]

On the other hand, from the formula (20) in Corollary A.1, we know that there exists a sub-sequence $\{n_{i}\uparrow\infty\}$ satisfying

[TABLE]

So when $n_{i}$ is big enough, we have

[TABLE]

Now we define a $p\times 1$ vector $\tilde{\beta}_{n_{i}}$ which includes only one non-zero component in the first coordinate, for $j=2,3,\cdots,p$

[TABLE]

where ${C_{0}}$ is a constant which will be specified later.

Then, when $n_{i}$ is large enough,

[TABLE]

Using the definition of $\tilde{\beta}_{n_{i}}$, condition (H.1)(ii) and the (42), we have

[TABLE]

Where $C_{1}$ is the determined constant. For ${R_{3{n_{i}}}}({{\tilde{\beta}}_{{n_{i}}}})$, we can see that

[TABLE]

By using Corollary A.2, which implies

[TABLE]

We also have sub-sequence $\{{{n^{\prime}}_{i}}{\rm{\}}}$ of $\{{{n}_{i}}{\rm{\}}}$ such that

[TABLE]

Analogous to the treatment of $R_{2n_{i}}(\tilde{\beta}_{n_{i}})$, we have

[TABLE]

Then we find that the sum of ${R_{2{{n^{\prime}}_{i}}}}({{\tilde{\beta}}_{{{n^{\prime}}_{i}}}})$ and ${R_{3{{n^{\prime}}_{i}}}}({{\tilde{\beta}}_{{{n^{\prime}}_{i}}}})$ is

[TABLE]

According to the $O(n)$ in (37), set $O({{n^{\prime}}_{i}})\leq C_{2}{{n^{\prime}}_{i}}$, and $C_{2}$ is a constant, we have

[TABLE]

Where $C_{3}$ is a constant. Therefore, according to (50) and (51), when ${n^{\prime}}_{i}$ is sufficiently large, we have

[TABLE]

Setting ${C_{3}}C_{0}^{2}-{C_{0}}(O(\sqrt{\frac{{\log\log{{n^{\prime}}_{i}}}}{{{{n^{\prime}}_{i}}}}})+C_{1})<0$, then ${C_{0}}<\frac{{O(\sqrt{\frac{{\log\log{{n^{\prime}}_{i}}}}{{{{n^{\prime}}_{i}}}}})+C_{1}}}{{{C_{3}}}}$. So we choose ${C_{0}}$ to be a certain positive constant which is smaller than $\frac{{O(\sqrt{\frac{{\log\log{{n^{\prime}}_{i}}}}{{{{n^{\prime}}_{i}}}}})+C_{1}}}{{{C_{3}}}}$. Thus,

[TABLE]

Notice the definition of $\hat{\beta}$ ($\hat{\beta}$ minimizes ${K_{n}}(\hat{\beta})$), which ensures that ${K_{{n^{\prime}}_{i}}}(\hat{\beta},{\beta_{0}})\leq{K_{{n^{\prime}}_{i}}}({{\tilde{\beta}}_{{{n^{\prime}}_{i}}}},{\beta_{0}})$. However when $n_{i}$ is sufficiently large, by assumption (36), we get

[TABLE]

via (40). This is a contradiction, thus the proof of Theorem 3.1 is completed. Therefore, for any $\alpha$ sub-model in $\Gamma_{c}$, $\mathop{\lim\sup}\limits_{n\to\infty}\frac{{\|{\hat{\beta}(\alpha)-{\beta_{0}}(\alpha)}\|}}{{\sqrt{{n^{-1}}\log\log n}}}=d\neq 0$  a.s. is true.

Appendix E Proof of Theorem 3.2

According to the proof of Theorem 3.1, from the definition of ${K_{n}}(\beta,{\beta_{0}})$, in order to prove (10) for any $\alpha$ sub-model in $\Gamma_{c}$, it is needed to prove the following

[TABLE]

Because $\hat{\beta}$ is the MLE of $\beta$, then ${l_{n}}(\beta_{0}(\alpha))\leq{l_{n}}(\hat{\beta}(\alpha))$, that is,

[TABLE]

Notice Lemma A.1 and Theorem 3.1

[TABLE]

Then, according to Corollary A.1 and Theorem 3.1, we have

[TABLE]

For ${R_{3n}}(\hat{\beta}(\alpha))$, by LIL in Corollary A.2, we get almost surely

[TABLE]

Therefore

[TABLE]

Combining the estimates (56), (57) and (59), we obtain

[TABLE]

Therefore, Theorem 3.2 is proved.

Appendix F Proof of Theorem 3.3

First, let ${{\hat{\beta}}^{*}}(\alpha)$ be the $p_{\alpha}$-dimensional vector which is defined by augmenting ${{\hat{\beta}}}(\alpha)$ with $p-p_{\alpha}$ 0’s such that the sub-vector of ${{\hat{\beta}}^{*}}(\alpha)$ indexed by $\alpha$ matches ${{\hat{\beta}}}(\alpha)$. Then, it is easy to see that proving the (11) is tantamount to give the proof of the following argument: for any incorrect model $\alpha\in\Gamma_{w}$

[TABLE]

Next, we define $l_{2}$ ball with radius $a:=\frac{1}{2}\min\limits_{1\leq i\leq p_{\alpha_{0}}}|\beta_{0}(\alpha_{0})_{i}|$,

[TABLE]

where $\alpha_{0}$ is the minimum dimension of the true model belonging to $\Gamma_{c}$, and $\beta_{0}(\alpha_{0})_{i}$ is the $i$th component of ${K_{n}}(\hat{\beta},{\beta_{0}})$.

Obviously, $B_{0}$ is a tight set. By the definition of incorrect model $\alpha\in\Gamma_{w}$, we get $\|\hat{\beta}^{*}(\alpha)-\beta_{0}\|\geq\min\limits_{1\leq i\leq p_{\alpha_{0}}}|\beta_{0}(\alpha_{0})_{i}|$. Then we have $\hat{\beta}^{*}(\alpha)\notin B_{0}$. By applying Theorem 3.1, for the large $n$, MLE $\hat{\beta}$ is almost surely an interior point of $B_{0}$. Applying the convexity of ${K_{n}}(\beta,{\beta_{0}})$, we have

[TABLE]

where $\partial B_{0}=\{\beta:\|\beta-\beta_{0}\|=\frac{1}{2}\min\limits_{1\leq i\leq p_{\alpha_{0}}}|\beta_{0}(\alpha_{0})_{i}|=a\}$, which is the boundary of $B_{0}$.

To prove (60), it is required to prove the following result:

[TABLE]

Using the similar argument in (27), we have

[TABLE]

Then

[TABLE]

Similarly, using (29) one has

[TABLE]

Then by taking supreme, we get

[TABLE]

For, using $\eqref{R3S}$ we similarly have

[TABLE]

By taking supreme, it gives

[TABLE]

According to (63), (65) and (66), we obtain

[TABLE]

Then

[TABLE]

Therefore, when $n\rightarrow\infty$, $\inf\inf\limits_{\beta\in\partial B_{0}}n^{-1}{K_{n}}(\beta_{0},{\beta})>0$. Hence (62) and (60) is established. Then Theorem 3.3 is proved.

Appendix G Proof of Theorem 3.4

The proof is to divide $\alpha$ into two cases. That is: 1.$\alpha\in\Gamma_{c}$; 2.$\alpha\in\Gamma_{w}$.

Step1: For any correct model $\alpha\in\Gamma_{c}$, employing Theorem 3.2 we have

[TABLE]

Minimizing over $\alpha$ gives

[TABLE]

Since the right hand side above should not be a decreasing function as $n$ is large, it leads to $C(n,\hat{\beta}(\alpha))>O(\log\log n)$. The reason is that for any correct model $\alpha\in\Gamma_{c}$ we need following equality:

[TABLE]

Note that the definition of ${\alpha_{0}}$ is ${\alpha_{0}}:=\mathop{\arg\min}\limits_{\alpha\in A,\alpha\in{\Gamma_{c}}}{p_{\alpha}}$. Then, by the definition of ${\hat{\alpha}}$ and the fact that $C(n,\hat{\beta}(\alpha))$ is an increasing function of ${p_{\alpha}}$, we have

[TABLE]

Thus $\hat{\alpha}={\alpha_{0}}$ a.s..

Step2: For any wrong model $\alpha\in\Gamma_{w}$, employing Theorem 3.3 we get

[TABLE]

Since $\alpha$ is wrong, thus ${\alpha_{0}}\neq\hat{\alpha}$, and then $\alpha_{0}$ should not be obtained by minimizing ${S_{n}}(\alpha)$. Therefore ${S_{n}}(\alpha)$ will increase as $\alpha$ approximates $\alpha_{0}$. Observe that minimizing ${S_{n}}(\alpha)$ is equivalent to minimizing $C(n,\hat{\beta}(\alpha))+O(n)$. Consequently, if we assume that minimizing $C(n,\hat{\beta}(\alpha))+O(n)$ is the same as minimizing $C(n,\hat{\beta}(\alpha))$, then this assumption leads to $\hat{\alpha}={\alpha_{0}}$ a contradiction with ${\alpha_{0}}\neq\hat{\alpha}$. Hence, we must have $C(n,\hat{\beta}(\alpha))<O(n)$.

In summary, if $O(\log\log n)<C(n,\hat{\beta}(\alpha))<O(n)$, then the model selection criteria are strongly consistent. Finally, via the order of AIC penalty term, it is not strongly consistent. Similarly, both the order of BIC and SCC are $\log n$ and it follows that they are strongly consistent.

Appendix H Proof of Theorem 3.5

Before the proof, we pose two LIL results for $\rho$-mixing process and $m$-dependent responses.

Lemma H.1**.**

*(LIL of $\rho$-mixing process, Corollary 9.2.1 in Lin and Lu (1997)) Let $\{Z_{n},n\geq 1\}$ be a strictly stationary $\rho$-mixing sequence with $E(Z_{1})=0$, $E(Z_{1}^{2})<\infty$. Let $S_{n}=\sum\limits_{k=1}^{n}Z_{k}$. Assume that

(i) $\sigma_{n}^{2}:=E(S_{n}^{2})\to\infty$ as $(n\rightarrow\infty)$;

(ii) $\rho(n)=O({(\log n)^{-1-\varepsilon}})$ for some $\varepsilon>0$.

Then, we have*

[TABLE]

Lemma H.2**.**

(LIL of $m$-dependent random variables, Chen (1997)) Let $\{Z_{n},n\geq 1\}$ be a real stationary strongly mixing sequences with $E(Z_{1})=0$, $E(Z_{1}^{2})<\infty$. If $\{Z_{n},n\geq 1\}$ is $m$-dependent and let $S_{n}=\sum\limits_{k=1}^{n}Z_{k}$, we have

[TABLE]

where $\sigma^{2}:=EZ_{1}^{2}+2\sum\limits_{k=2}^{m+1}{E{X_{1}}{X_{k}}}$.

The proof of (15) for $\rho$-mixing responses or $m$-dependent responses is slightly different from what we mentioned in Theorem 3.1 concerning independence case, and the main dissimilarity is the calculating of the variance of the score function. We cannot ignore the non-zero cross terms which will be included as the covariance parts.

(i): $\rho$-mixing sequence: For $j=1,\cdots,p$, let

[TABLE]

where {{\mathord{\buildrel{\lower 3.0pt\hbox{\scriptscriptstyle\smile}}\over{w}}}_{k}}:={w_{k}}{x_{kj}}\dot{u}(x_{k}^{t}\beta_{0}).

We have

[TABLE]

Lemma H.3**.**

(Davydov’s inequality) Let $X\in\mathcal{F}_{1}^{k},Y\in\mathcal{F}_{k+n}^{\infty}$ with $E|X|^{p}<\infty$ and $E|Y|^{p}<\infty$ (${p^{-1}}+{q^{-1}}<1$), then

[TABLE]

(see Lemma 1.2.4 in Lin and Lu (1997) or Corollary 1.1.1 in Bosq (1998)).

By Davydov’s inequality above, let $p=q=3$ in Lemma H.3, we have

[TABLE]

By (H.5): The error sequence ${\varepsilon_{i}}$ satisfies $\rho$-mixing condition with geometric decay: $\rho(m)=O({r^{-m}})$, and the relationship between $\alpha$-mixing coefficient and $\rho$-mixing coefficient meets: $\alpha(n)\leq\frac{{\rm{1}}}{{\rm{4}}}\rho(n)$, we have

[TABLE]

Then (69) and (70) imply ${I_{n}}(\beta)(k,k)=O(n)\to\infty$. So condition (i) in Lemma H.1 is satisfied. And condition (ii) in Lemma H.1 is also valid since the (H.5) gives $\rho(n)=O({r^{-n}})=O({(\log n)^{-1-\varepsilon}})$.

(ii) $m$-dependent sequence: For $j=1,\cdots,p$, we have

[TABLE]

Under (H.6), then we can use a result of the LIL for real stationary strongly mixing sequences with $m$-dependent properties, see Lemma H.2.

In summary, the LIL for the $\alpha$-mixing or $m$-dependent score functions is

[TABLE]

Therefore, the proof of Theorem 3.5 can be imitated from the proof of Theorem 3.1.

After Observing the proofs of Theorem 3.2, Theorem 3.3 and Theorem 3.4 do not involve the calculation of the second moment of the weighted score function, we conclude that the corresponding new proofs are consistent with the proofs of the case of independent responses.

We only pay attention to the equations which concern the LIL of the several weighted score functions, i.e. (31), (41) and (47). Next, for weakly dependent version of Theorem 3.3, the LIL is also true when the weighted score function is weakly dependent sums. As for Theorem 3.3, the proof is the same by applying the LIL of weakly dependent weighted score function.

Then as for two dependent cases, assuming that conditions (H.1) to (H.6) are satisfied, we also have the same conclusions as in Theorem 3.2 and Theorem 3.3. The proofs in other places in Theorem 3.6 are consistent with Theorem 3.4. Therefore we do not make repetitions here.