5_xxx_qian.dvi EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS Vol. 3, No. 3, 2010, 417-434 ISSN 1307-5543 – www.ejpam.com SPECIAL ISSUE ON GRANGER ECONOMETRICS AND STATISTICAL MODELING DEDICATED TO THE MEMORY OF PROF. SIR CLIVE W.J. GRANGER Law of Iterated Logarithm and Strong Consistency in Poisson Regression Model Selection Guogi Qian Department of Mathematics and Statistics, University of Melbourne, VIC 3010, Australia Abstract. In this paper we first derive a law of iterated logarithm for the maximum likelihood esti- mator of the parameters in a Poisson regression model. We then use this result to establish the strong consistency of a class of model selection criteria in Poisson regression model selection. We show that under some general conditions, a model selection criterion, which consists of a minus maximum log- likelihood and a penalty term, will select the simplest correct model almost surely if the penalty term increases with model dimension and has an order in between O(log log n) and O(n). 2000 Mathematics Subject Classifications: 62F12, 62J12, 60F15 Key Words and Phrases: Law of iterated logarithm; Poisson regression; Maximum likelihood estima- tor; Model selection; Strong consistency 1. Introduction Poisson regression model is a widely used econometric and statistical tool for studying the relationship between a Poisson-type response variable and a set of explanatory variables. A familiar example is the analysis of contingency tables of categorical data. In addition to parameter estimation, another important inference task in Poisson regression is searching for a subset of available explanatory variables that can best explain or predict the response. This amounts to the Poisson regression model or variable selection. Many papers can be found in recent literature in the area of model selection, which deal with different models in different ways. We refer to [4] and [14] and references therein for the detailed survey. It appears that, while it might be implied according to some general model selection prin- ciple, the Poisson regression model selection method by itself has hardly been investigated in Email address: g.qian�ms.unimelb.edu.au http://www.ejpam.com 417 c© 2010 EJPAM All rights reserved. G. Qian / Eur. J. Pure Appl. Math, 3 (2010), 417-434 418 a formal and rigorous way. The lack of a formal theory for Poisson regression model selec- tion creates uncertainty and inconvenience for people applying the method for practice. This motivates the writing of this paper which focuses on the asymptotic performance of a class of model selection criteria including AIC, BIC, Mallows Cp and stochastic complexity or min- imum description length for Poisson regression models. A byproduct of this asymptotic study is the establishment of the law of iterated logarithm for the maximum likelihood estimators (MLE) in the Poisson regression models. The convergence rate of the MLE provided by the law of iterated logarithm is very useful in deriving precise approximations for likelihood function. In the paper we first set up a model selection framework for Poisson regression models and review several general model selection criteria such as AIC [1], BIC [18] and stochastic complexity or minimum description length [16] in section 2. In section 3 we present the main results and the conditions for ensuring these results. We have shown that when the employed model is a correct one, the MLE β̂ converges almost surely to the true parameter value β0 with a rate not slower than O( p n−1 log log n). We have also shown that, for a model selection criterion consisting of the minus log-likelihood and a penalty term, it will select the simplest correct model almost surely if the penalty term is an increasing function of the model dimension and is of an order higher than O(log log n) and lower than O(n). The detailed proof of these results are given in section 4 and the appendix. The paper is concluded with a discussion given in section 5. 2. Model Selection in Poisson Regression Model The problem to our interest is whether any component of a given explanatory vector x = (x1, · · · , xp) t has any effect on a response variable Y . When Y is a count variable, it is often sensible to assume a Poisson distribution for Y which has a probability function P(Y = y) = (y!)−1µy e−µ (y = 0,1,2, · · · ). Then the problem can be studied in the framework of a log-linear regression model which assumes a linear predictor η = xtβ for logarithm of the mean of Y , i.e., logµ = η = xtβ , where β = (β1, · · · ,βp) t is the unknown parameter vector of interest. Now let Yn = (y1, · · · , yn) t be the n independent observations from Y , with the corre- sponding explanatory vectors being x1, · · · ,xn. Denote Xn = (x1, · · · ,xn) t as the design matrix. Then under the log-linear regression model considered, the distribution for yi is Poisson(µi) with µi = eηi = ext i β ; and the log-likelihood function for the parameter β is ℓ(β |Yn, Xn) = n ∑ i=1 {− log yi!+ yi logµi −µi}= − n ∑ i=1 log yi!+ n ∑ i=1 {yix t iβ − ext i β}. (1) The maximum likelihood estimator(MLE) β̂ is defined to be β̂ = arg max β ℓ(β |Yn, Xn) = arg min β n ∑ i=1 {ext i β − yix t iβ} G. Qian / Eur. J. Pure Appl. Math, 3 (2010), 417-434 419 which can be solved from ∂ ℓ ∂ β = n ∑ i=1 (yi −µi)xi = 0. In this paper we will show that the estimation error ||β̂ − β0|| = O( p n−1 log log n) almost surely under some general conditions. Here β0 is the true value of β and || · || is the Euclidean norm. The Fisher information for the parameter β can be found to be In(β) = −E ∂ 2ℓ ∂ β∂ β t = − ∂ 2ℓ ∂ β∂ β t = X t nU nXn, (2) where Un = diag{µ1, · · · ,µn}. With the Poisson regression model logµ = xtβ , the effect of each x variable on Y can be measured by the value of the corresponding β component. Thus if any of the β compo- nents equals 0 or is close to 0, there would be no necessity to include the corresponding x components into the model. Since the values of β can only be estimated, one needs a statisti- cal model selection criterion to determine which x components have significant effects on Y thus should be included in the model. The maximum likelihood principle cannot serve as a model selection criterion because the maximum likelihood for the full model including all the available x components is always greater than the maximum likelihood for a sub-model using a subset of the x components. But a model selection criterion can be based on a penalised log-likelihood. Let α be a pα-component sub-vector of (1,2, · · · , p). Let xα and β(α) be the sub-vectors of x and β indexed by α respectively. Further let logµα = ηα = xt αβα be a Poisson regression model containing a subset of explanatory variables given by xα. The penalised log-likelihood based model selection criterion can be expressed as S(ηα) = −ℓ(β̂(α)|Yn, Xnα) + C(n, β̂(α)), (3) where the first term is the minus maximum log-likelihood measuring the goodness of fit of model ηα, and the second term measures the complexity of the model. The matrix Xnα comprises those columns of Xn indexed by α; and β̂(α) is the MLE of β(α). Under the criterion (3), those sub-models having both better goodness of fit and smaller complexity will be preferred than the others; and the best model will be the one achieving the small- est S(ηα) value. Many commonly used model selection criteria, such as AIC [1], BIC [18], Cp [7] and stochastic complexity criterion (SCC)[16, 17, 11], are of the form given by (3). For example, for AIC and Cp C(n, β̂(α)) = pα; for BIC C(n, β̂(α)) = 1 2 pα log n; and for SCC C(n, β̂(α)) = 1 2 log |In(β̂(α))|+ ∑pα i=2 log(|β̂(α)i|+ ǫn−1/4) where β̂(α)i is the i-th component of β̂(α), and ǫ is a specified quantity to ensure the invariance of the SCC [see 11, for details]. Assuming that the model logµ = xtβ is the full model which includes all the explanatory variables available and the first component of x is an intercept term, there will be in total 2p−1 sub-models of the form logµα = xt αβ(α) for selection, provided that only those models having an intercept term are considered. In this paper we assume that some components of β0 = (β01, · · · ,β0p) t , the true value of β , are equal to 0. We also use α to represent a Poisson regression sub-model logµα = xt αβ(α), which is actually a one-to-one representation. Then all the 2p − 1 sub-models can be classified into the following two groups: G. Qian / Eur. J. Pure Appl. Math, 3 (2010), 417-434 420 1. Ac = {α : β0i = 0 for any i 6∈ α}; 2. Aw = {α : β0i 6= 0 for some i 6∈ α}. Apparently, every sub-model in Aw is a wrong model which misses at least one x variable having non-zero effect on Y , while every sub-model inAc is a correct model which includes all x variables having non-zero effects on Y . But the models inAc may contain some redundant x variables having no effects on Y . An ideal model selection criterion should render the selection of the simplest correct model in Ac containing no redundant explanatory variables. In this paper we will show that the penalised log-likelihood based model selection criterion, under some general conditions, selects the simplest correct model with probability 1 as the sample size n goes to infinity. For simplicity of the presentation, we assume the simplest correct model in Ac to be unique, which is the case if all components of x are linearly independent of each other. 3. Conditions and Main Results Let λ1{S} ≤ · · · ≤ λp{S} be the p eigenvalues of a p × p symmetric matrix S. Also let b = 1 2 min1≤i≤pα0 |β0(α0)i|, where α0 is the correct model in Ac with the smallest dimension, and β0(α0)i is the i-th component of β0(α0). We assume b > 0 in this paper. Note that b is only used in the proof of Theorem 3 in this paper; and b = 0 represents the case whereAw is an empty set thereby Theorem 3 is not applicable. Further define δn = q max 1≤i≤n µ0ix t i In(β0) −1xi and ξn = q max 1≤i≤n xt i In(β0) −1xi where µ0i = ext i β0 is the true value of µi. The following conditions will be required in various places in proving our main results: (C.1). limn→∞λ j{In(β0)} = ∞, j = 1, · · · , p. Also there exists a constant b0 such that 0 < λp{In(β0)} ≤ b0λ1{In(β0)}. (C.2). b1n≤ λp{In(β0)} ≤ b2n for some positive constants b1 and b2. (C.3). ξn p log logλp{In(β0)}= o(1). (C.4). δn p logλp{In(β0)} p log logλp{In(β0)} = o(1). (C.5). δn logλp{In(β0)} p log logλp{In(β0)} = o(1). (C.6). ξn logλp{In(β0)} p log logλp{In(β0)} = o(1). (C.7). λ1{X t nMnXn} ≥ b3n for some constant b3 > 0, where Mn = diag{µ01e−2b‖x1‖, · · · ,µ0ne−2b‖xn‖}. G. Qian / Eur. J. Pure Appl. Math, 3 (2010), 417-434 421 Note that these conditions are not completely independent of each other. Firstly, one can see that both (C.1) and (C.2) are implied from b1n ≤ λ1{In(β0)} ≤ λp{In(β0)} ≤ b2n. Sec- ondly, conditions (C.1), (C.2) and (C.7) together suggest that all the eigenvalues of In(β0) and XnMnXn are of order O(n). Thirdly, conditions (C.3) to (C.6) together with (C.2) indicate ξn p log log n→ 0, δn p log n log log n→ 0, δn log n p log log n→ 0 and ξn log n p log log n→ 0 respectively. Finally, under condition (C.1), (C.3) is implied by (C.6); and (C.4) is implied by (C.5). Although conditions (C.1) to (C.7) may be simplified according to the preceding discussion, we prefer not to do so in order to clarify that to what extent each condition is required in the proof. The conditions (C.1) to (C.7) are essentially about the behavior of the explanatory vari- ables x. Roughly speaking, they mean most of the observations {x1, · · · ,xn} should be finite and stay away from 0; and if a subsequence of {x1, · · · ,xn} diverges to infinity, it should do so with an appropriate rate. In fact, if we assume x is a random vector and x1, · · · ,xn as i.i.d. observations from x, then the following are sufficient for (C.1) to (C.7) to hold: (S.1). The moment generating function Eext s exists for ||s|| ≤ ||β0|| + s0 for some constant s0 > 0. This implies that all of E(extβ0xtx)κ, E(extβ0−2b‖x‖xtx)κ and E(xtx)κ are finite for some κ > 1. (S.2). P(xtv 6= 0) > 0 for all v 6= 0 in R p, which implies Eextβ0xxt , Eextβ0−2b‖x‖xxt and Exxt are all positive definite. To see the sufficiency of (S.1) and (S.2), one can apply the strong law of large numbers for the i.i.d. random variables x1, · · · ,xn, · · · under condition (S.1), which gives the following results: 1 n X t nUnXn− Eextβ0xxt a.s.→ 0, 1 n X t nMnXn− Eextβ0−2b‖x‖xxt a.s.→ 0. These results together with (S.2) imply (C.1),(C.2) and (C.7). The conditions (C.3) to (C.6) are implied from (C.1), (C.2) and the fact that, under (S.1) δ2(1+κ′) n ≤ λ1{In(β0)}−1−κ′ n ∑ j=1 e xt j β0(1+κ ′) (xt jx j) 1+κ′ = O(n−κ ′ ) a.s. and ξ2(1+κ′) n ≤ λ1{In(β0)}−1−κ′ n ∑ j=1 (xt jx j) 1+κ′ = O(n−κ ′ ) a.s. for some κ′ > 0. In this paper we will regard the observations x1, · · · ,xn as deterministic for simplicity of the presentation. There is no essential complication with random xi ’s. In this paper we have obtained the following results. Theorem 1. Suppose conditions (C .1) to (C .6) are satisfied. Then for any correct model α ∈ Ac, ||β̂(α)− β0(α)||= O( p n−1 log log n) a.s.. (4) G. Qian / Eur. J. Pure Appl. Math, 3 (2010), 417-434 422 Further, there exists a constant d > 0 such that for any α ∈Ac lim sup n→∞ ||β̂(α)−β0(α)|| p n−1 log log n = d a.s.. (5) Hence the M LE β̂(α) follows the law of iterated logarithm. Theorem 2. Under conditions (C .1) to (C .6), for any correct model α ∈Ac , 0≤ ℓ(β̂(α)|Yn, Xnα)− ℓ(β0|Yn, Xn) = O(log log n) a.s.. (6) Theorem 3. Under conditions (C .1) to (C .7), for any incorrect model α ∈ Aw, we have lim sup n→∞ 1 n {ℓ(β̂(α)|Yn, Xnα)− ℓ(β0|Yn, Xn)} < 0 a.s.. (7) From Theorems 2 and 3 we know that the maximum log-likelihood of any correct model is almost surely greater than the unknown true log-likelihood of the full model by an amount bounded by |O(log log n)|. On the other hand, the maximum log-likelihood of any incorrect model in Aw is almost surely smaller than the true log-likelihood of the full model by an amount greater than τn with τ > 0 when n is sufficiently large. Therefore, if we use a penalised log-likelihood based criterion of form (3) for model selection, we will almost surely select the simplest correct model inAc if the penalty term C(n, β̂(α)) is an increasing function of the model dimension pα and is of an order in between O(log log n) and O(n). We call a model selection criterion strongly consistent if it selects the simplest correct model almost surely; and consistent if almost surely it only selects one of the correct models. From this discussion we have the following: Theorem 4. For a Poisson regression model satisfying conditions (C .1) to (C .7), both model selection criteria BIC and SCC are strongly consistent, while AIC is consistent but not necessarily strongly consistent. Proof. As the Fisher information’s determinant |I(β(α))| is typically of order O(npα), both the penalty terms of SCC and BIC are increasing functions of the model dimension pα and are of order O(log n), it follows from Theorem 2 and Theorem 3 that both SCC and BIC are strongly consistent. AIC is not necessarily strongly consistent because its penalty term is of order O(1). But AIC is clearly consistent because its criterion value for a correct model is almost surely smaller than that for any incorrect model by an amount greater than τn when n is sufficiently large. The proof of Theorems 1 to 3 will be the focus of the next section. 4. Proof of the Results The key to proving our main results lies on the convexity and quadratic approximation of the negative log-likelihood function, the Normal and Gamma approximations of the Poisson G. Qian / Eur. J. Pure Appl. Math, 3 (2010), 417-434 423 probabilities, and the law of iterated logarithm for independent random variables. The idea of using the convexity property is broadly seen in establishing asymptotic representations of the M-estimators in linear models, see e.g. [15, 20, 12] among the others. By the definition of ξn and conditions (C.2) and (C.3) it is easy to see that there exists a sequence of positive numbers {τn} satisfying: τn ↑ ∞, τnξn p log log n→ 0 and τn p n−1 log log n ↓ 0. Using τn we introduce two sequences of subsets: An = {β : ||β − β0|| ≤ τn p n−1 log log n} ∂ An = {β : ||β − β0||= τn p n−1 log log n}. It is clear that A1 ⊃ A2 ⊃ A3 ⊃ · · · ⊃ An. Further we define H(β , n) = ℓ(β0|Yn, Xn)− ℓ(β |Yn, Xn) = n ∑ k=1 {ext k β − ext k β0 − ykxt k(β − β0)}, and K(t, s) = et − es − es(t − s). By these definitions it follows that H(β , n) = n ∑ k=1 K(xt kβ ,xt kβ0)− n ∑ k=1 (yk −µ0k)x t k(β − β0) def = R1(β , n) + R2(β , n). (8) Before proving the main results we need to establish some preliminary results. Lemma 1. The function K(t, s) defined has the following properties: (i). K(t, s) ≥ 0 for any real numbers t and s. (ii). K(t, s) is strictly convex with respect to t. (iii). For any ∆ > 0, 1 2 es−2∆(t − s)2 ≤ K(t, s) ≤ 1 2 es+2∆(t − s)2 if |t − s| ≤∆. The proof of Lemma 1 will be give in Appendix. Lemma 2. Let W be a Poisson(θ) random variable. Then for any w ≥ 0 the following inequalities hold: P(W ≤ w)≤ (2π)− 1 2 ∫ (w+1−θ )/pθ −∞ e− 1 2 t2 d t, (9) P(W ≤ w) ≥ [Γ(θ + 1)]−1 ∫ w 0 tθ e−t d t. (10) G. Qian / Eur. J. Pure Appl. Math, 3 (2010), 417-434 424 The results of Lemma 2 were obtained by [2] which can also be found in [6, p. 102]. Lemma 3 (Law of the Iterated Logarithm). Let {Zn, n ≥ 1} be independent random vari- ables with EZn = 0, EZ2 n = σ 2 n and s2 n = ∑n k=1σ 2 k → ∞. If |Zn| ≤ dn a.s., where dn = o((s2 n/ log log s2 n) 1/2), then lim sup n→∞ ±∑n k=1 Zk p 2s2 n log log s2 n = 1 a.s.. This lemma and its proof can be found in e.g. [3, pp. 373-374] and [9, pp. 239-246]. Lemma 4. Under conditions (C.1), (C.2) and (C.4) to (C.6), we have lim sup n→∞ ±∑n k=1(yk −µ0k)xk j p 2In(β0)( j, j) log log In(β0)( j, j) = 1 a.s. for j = 1, · · · , p. (11) Here xk j is the j-th element of xk and In(β0)( j, j) is the ( j, j)-th element of In(β0). Equation (11) suggests that {(yk−µ0k)xk j, k = 1,2, · · · } obeys the law of iterated logarithm. Accordingly, we have ∂ ℓ ∂ β |β=β0 = n ∑ k=1 (yk −µ0k)xk = X t n(Yn−µ0) = O( p n log log n) a.s. (12) where µ0 = (µ01, · · · ,µ0n) t is the true mean vector. Proof. The result (12) is obvious from (11) and condition (C.2). Hence we only need to prove (11). Without losing generality we assume all xk j > 0. Using the information that yk ∼ Poisson(µ0k) and the definition of In(β0) it is easy to verify that for j = 1, · · · , p E(yk −µ0k)xk j = 0, (13) n ∑ k=1 E((yk −µ0k)xk j) 2 = n ∑ k=1 µ0k x2 k j = In(β0)( j, j)→∞ (14) as n→∞ by condition (C.1). From now on we proceed to show that for j = 1, · · · , p |(yn−µ0n)xnj| ≤ o(dnj) a.s. (15) where dnj = p In(β0)( j, j)/ log log In(β0)( j, j)→ ∞ by (14). For any ǫ > 0, it is easy to see that P{|(yn −µ0n)xnj| > ǫdnj} ≤ P{yn > µ0n+ ǫdnj x −1 nj }+ P{yn < µ0n − ǫdnj x −1 nj }. (16) Applying (9) of Lemma 2 we have P{yn < µ0n − ǫdnj x −1 nj } ≤ (2π)− 1 2 ∫ (1−ǫdn j x −1 n j )/ p µ0n −∞ e− 1 2 t2 d t. (17) G. Qian / Eur. J. Pure Appl. Math, 3 (2010), 417-434 425 Note that from condition (C.1) and the inequality λ1{In(β0)} ≤ In(β0)( j, j) ≤ λp{In(β0)} we have µ0n x2 nj ≤ µ0nxt nxn ≤ λp{In(β0)}µ0nxt n In(β0) −1xn ≤ λp{In(β0)}δ2 n ≤ λp{In(β0)} λ1{In(β0)} In(β0)( j, j) log log In(β0)( j, j) δ2 n log logλp{In(β0)} ≤ b0d2 njδ 2 n log logλp{In(β0)}. (18) Thus when µ0n ≥ 1, 1− ǫdnj x −1 njp µ0n = 1 p µ0n − ǫ √ √ √ √ d2 nj µ0n x2 nj ≤ 1− ǫ p b0δ 2 n log logλp{In(β0)} → −∞ by condition (C.4). This implies that 1− ǫdnj x −1 njp µ0n ≤ −1 2 ǫ(b0δ 2 n log logλp{In(β0)})− 1 2 when n is sufficiently large. (19) From (17), (19) and a well-known inequality ∫ ∞ a e− 1 2 t2 d t < 1 a e− 1 2 a2 for all a > 0 [see 3, p.49], it follows that when µ0n ≥ 1 and n is sufficiently large, P{yn < µ0n − ǫdnj x −1 nj } ≤ (2π)− 1 2 ∫ − 1 2 ǫ(b0δ 2 n log logλp{In(β0)})− 1 2 −∞ e− 1 2 t2 d t < (2π)−12ǫ−1(b0δ 2 n log logλp{In(β0)}) 1 2 e− 1 8 ǫ2(b0δ 2 n log logλp{In(β0)})−1 ≤ (logλp{In(β0)})− 1 2 (λp{In(β0)})−2 (20) where the last inequality follows from condition (C.4). When µ0n < 1, it follows from (18) and conditions (C.1) and (C.4) that µ0n − ǫdnj x −1 nj = p µ0n( p µ0n− ǫ √ √ √ √ d2 nj µ0n x2 nj )≤pµ0n(1− ǫ(b0δ 2 n log logλp{In(β0)})− 1 2 )< 0 if n is sufficiently large, which suggests P{yn < µ0n − ǫdnj x −1 nj } = 0 if n is sufficiently large. Therefore the result (20) is true for any µ0n > 0, which implies that ∞ ∑ n=1 P{yn < µ0n − ǫdnj x −1 nj }<∞ (21) G. Qian / Eur. J. Pure Appl. Math, 3 (2010), 417-434 426 if further condition (C.2) holds. Now let {cn : 0 < cn ≤ 1 2 , n = 1,2, · · · } be a sequence of numbers to be determined. Applying the result (10) of Lemma 2 and the property for the Gamma function [Γ(θ + 1)]−1 ∫∞ 0 tθ e−t d t = 1, we have P{yn > µ0n+ ǫdnj x −1 nj } ≤ [Γ(µ0n+ 1)]−1 ∫ ∞ µ0n+ǫdn j x −1 n j tµ0n e−t d t ≤ e −cn(µ0n+ǫdn j x −1 n j ) [Γ(µ0n+ 1)]−1 ∫ ∞ µ0n+ǫdn j x −1 n j tµ0n e−(1−cn)t d t ≤ e −cn(µ0n+ǫdn j x −1 n j ) (1− cn) −(µ0n+1) = 1 1− cn � e−cn 1− cn �µ0n e −ǫcndn j x −1 n j ≤ 2 1− cn + 1 2 c2 n 1− cn !µ0n e −ǫcndn j x −1 n j = 2 � 1+ c2 n 2(1− cn) �µ0n e −ǫcndn j x −1 n j . (22) We take cn = min{1 2 ,µ − 1 2 0n } so c2 n 2(1−cn) ≤ min{1 4 ,µ−1 0n }. By considering the two cases µ0n ≤ 4 and µ0n > 4 separately and using the property that (1+ 1 a )a ↑ e as a ↑ ∞, it is easy to see that � 1+ c2 n 2(1− cn) �µ0n ≤max{ � 5 4 �4 , e} = e. (23) Applying (18) one can show that when n is sufficiently large ǫdnjµ − 1 2 0n x−1 nj ≥ ǫb − 1 2 0 (δn p log logλp{In(β0)})−1 ≥ 2 logλp{In(β0)} (24) under conditions (C.1) and (C.5). In the same way as proving (18) one can show that under condition (C.1) x2 nj ≤ b0d2 njξ 2 n log logλp{In(β0)}. (25) By (25) and condition (C.6) it follows that 1 2 ǫdnj x −1 nj ≥ 1 2 ǫb − 1 2 0 (ξn p log logλp{In(β0)})−1 ≥ 2 logλp{In(β0)} (26) when n is sufficiently large. By (24), (26) and the fact that cn = min{1 2 ,µ − 1 2 0n } it follows that when n is sufficiently large, ǫcndnj x −1 nj ≥ 2 logλp{In(β0)}. (27) Now from (22), (23) and (27) we have P{yn > µ0n+ ǫdnj x −1 nj } ≤ 2eλp{In(β0)}−2 (28) G. Qian / Eur. J. Pure Appl. Math, 3 (2010), 417-434 427 when n is sufficiently large. From (28) and condition (C.2) it follows that ∞ ∑ n=1 P{yn > µ0n+ ǫdnj x −1 nj }<∞. (29) Following the results (16), (21) and (29) we have ∞ ∑ n=1 P{|(yn −µ0n)xnj| > ǫdnj}<∞. Hence by the Borel-Cantelli lemma, P{|(yn−µ0n)xnj| > ǫdnj occurs infinitely often} = 0 for any ǫ > 0, which implies that (15) is true. Since (13) to (15) are true, the result (11) is followed by applying Lemma 3 for the independent random variables {(yk −µ0k)xk j, k = 1,2, · · · }. Proof. (Theorem 1) Clearly it is sufficient to prove (4) only for the full model: ||β̂ −β0||= O( p n−1 log log n) a.s. (30) Applying result (iii) of Lemma 1 with t = xt k β , s = xt k β0 and ∆= |xt k β − xt k β0|, it follows that K(xt k β ,xt k β0)≥ 1 2 e−2|xt k (β−β0)|µ0k[x t k (β −β0)] 2. (31) Following the definition of ξn and condition (C.2) one can find that max 1≤k≤n ||xk||2 ≤ λp{In(β0)} max 1≤k≤n xt k In(β0) −1xk = λp{In(β0)}ξ2 n ≤ b2nξ2 n. (32) Thus by (32) and Cauchy-Schwarz inequality, max 1≤k≤n |xt k (β−β0)|I(β∈∂An) ≤ max 1≤k≤n ||xk|| · ||β−β0||I(β∈∂An)≤ p b2ξnτn p log log n (33) where I(β ∈ ∂An) is an indicator function indicating that only those β in ∂An will be under consideration. (This type of definition for the indicator function will be used in the rest of the paper.) It follows from (8) and (31) to (33) that R1(β , n)I(β∈∂An)≥ 1 2 e−2max1≤k≤n |xt k (β−β0)| n ∑ k=1 µ0k[x t k(β − β0)] 2I(β∈∂An) ≥ 1 2 e−2 p b2ξnτn p log log n(β − β0) t In(β0)(β − β0)I(β∈∂An) ≥ 1 2 e−2 p b2ξnτn p log log nλ1{In(β0)}||β − β0||2I(β∈∂An). (34) G. Qian / Eur. J. Pure Appl. Math, 3 (2010), 417-434 428 By conditions (C.1) and (C.2) and the fact that τnξn p log log n→ 0, it comes after (34) that there exists a constant b4 > 0 such that R1(β , n)I(β∈∂An)≥ b4τ 2 n log log n. (35) On the other hand, by result (12) of Lemma 4 and (8), |R2(β , n)|I(β∈∂An)≤ || n ∑ k=1 (yk −µ0k)xk|| · ||β − β0||I(β∈∂An) = O( p n log log n)τn p n−1 log log n= O(1)τn log log n a.s.. (36) Knowing (35) and (36) we can find a constant b5 > 0 so that H(β , n)I(β∈∂An) = {R1(β , n) + R2(β , n)}I(β∈∂An)≥ b5τ 2 n log log n a.s.. (37) It is easy to see that H(β , n) is convex by Lemma 1 and that H(β0, n) = 0. This and (37) suggests that the MLE β̂ which also minimizes H(β , n) must be inside the subset An almost surely; namely ||β̂ − β0|| ≤ τn p n−1 log log n a.s.. (38) Equation (38) implies (30) because the sequence {τn} can be chosen to diverge as slowly as possible. We now proceed to prove (5) for the full model. Suppose (5) does not hold for the full model. This implies that ||β̂ − β0|| = o( p n−1 log log n) a.s., (39) knowing that (30) is true. By applying result (iii) of Lemma 1, (8) and (32) it follows that R1(β̂ , n)≤ 1 2 e2max1≤k≤n |xt k (β̂−β0)| n ∑ k=1 µ0k[x t k(β̂ − β0)] 2 ≤ 1 2 e2 p b2nξn||β̂−β0||(β̂−β0) t In(β0)(β̂−β0)≤ 1 2 e2 p b2nξn||β̂−β0||λp{In(β0)}||β̂−β0||2. (40) Thus by (39) and conditions (C.2) and (C.3) we have R1(β̂ , n) = o(1) log log n a.s.. Corre- sponding to (36) it can be seen that R2(β̂ , n) = o(1) log log n a.s. under the assumption of (39). Hence we have that under assumption (39) H(β̂ , n) = R1(β̂ , n) +R2(β̂ , n) = o(1) log log n a.s.. (41) On the other hand, from (11) of Lemma 4 we know there exists a sequence of positive integers {ni ↑ ∞} such that lim i→∞ ∑ni k=1 (yk −µ0k)xk1 p 2Ini (β0)(1,1) log log Ini (β0)(1,1) = 1 a.s.. G. Qian / Eur. J. Pure Appl. Math, 3 (2010), 417-434 429 Thus when ni is sufficiently large, ∑ni k=1 (yk −µ0k)xk1 p 2Ini (β0)(1,1) log log Ini (β0)(1,1) ≥ 1 2 a.s.. (42) Now define a p× 1 vector β̃ni with β̃ni ( j) = β0 j for j = 2, · · · , p and β̃ni (1) = b1 4b0 b2 È 2 log log Ini (β0)(1,1) Ini (β0)(1,1) + β01. Then from (42) it follows that, when ni is sufficiently large R2(β̃ni , ni) = ni ∑ k=1 (µ0k− yk)x t k (β̃ni −β0) = ni ∑ k=1 (µ0k− yk)xk1(β̃ni (1)− β01) ≤ −1 2 p 2Ini (β0)(1,1) log log Ini (β0)(1,1) · b1 4b0 b2 È 2 log log Ini (β0)(1,1) Ini (β0)(1,1) = − b1 4b0 b2 log log Ini (β0)(1,1) a.s.. (43) Note that by conditions (C.1) and (C.2) b2n≥ λp{In(β0)} ≥ In(β0)(1,1)≥ λ1{In(β0)} ≥ b1 b0 n (44) and accordingly 2 log log n≥ log log In(β0)(1,1)≥ 1 2 log log n when n is sufficiently large. (45) It follows that p b1 2 p b0 b2 Æ n−1 i log log ni ≥ β̃ni (1)− β01 ≥ b1 4b0 b2 p b2 Æ n−1 i log log ni (46) when ni is sufficiently large. Using (46), (32) and the fact that ξn p log log n → 0 one can show that max 1≤k≤ni |xt k (β̃ni−β0)| ≤ max 1≤k≤ni ||xk||(β̃ni (1)−β01) ≤ 1 2 È b1 b0 b2 ξni p log log ni ≤ 1 2 log2 (47) when ni is sufficiently large. Similar to proving (40), by (44) and (47) one can see that R1(β̃ni , ni)≤ λp{Ini (β0)}||β̃ni − β0||2 ≤ b2ni b2 1 16b2 0 b2 2 2 log log Ini (β0)(1,1) Ini (β0)(1,1) G. Qian / Eur. J. Pure Appl. Math, 3 (2010), 417-434 430 ≤ b1 8b0 b2 log log Ini (β0)(1,1) when ni is sufficiently large. (48) Thus, by (43), (48) and (45) it follows that when ni is sufficiently large, H(β̃ni , ni) = R1(β̃ni , ni) + R2(β̃ni , ni) ≤ − b1 8b0b2 log log Ini (β0)(1,1)≤ − b1 16b0b2 log log ni a.s.. (49) Since β̂ ≡ β̂n is the MLE that minimizes H(β , n), inferring from (49) we have H(β̂ni , ni)≤ H(β̃ni , ni)≤ − b1 16b0b2 log log ni a.s. (50) when ni is sufficiently large. Clearly, (50) is contradictory to (41) which suggests that (39) is wrong. Therefore (5) is true for the full model and consequently so for the other correct models inAc. Proof. (Theorem 2) As in proving Theorem 1, we only need to prove (6) for the full model which is equivalent to 0≥ H(β̂ , n) = R1(β̂ , n) + R2(β̂ , n) = O(log log n) a.s. (51) by the definition of H(β , n). The inequality part of (51) is obvious because β̂ is the MLE of β . Note that result (40) is also valid here, so by Theorem 1 and conditions (C.2) and (C.3) 0≤ R1(β̂ , n)≤ 1 2 e2 p b2nξn||β̂−β0||λp{In(β0)}||β̂−β0||2 = O(log log n) a.s.. (52) On the other hand, by result (12) of Lemma 4 and (4) of Theorem 1 we have |R2(β̂ , n)| ≤ || n ∑ k=1 (yk −µ0k)xk|| · ||β̂ − β0|| = O(log log n) a.s.. (53) By (52) and (53) it follows that |H(β̂ , n)| = O(log log n) a.s. which suffices the proof of the theorem. Proof. (Theorem 3) First we extend the pα×1 vector β̂(α), the MLE of βα, to a p×1 vector β̂∗(α) by inserting p− pα 0’s into β̂(α) in such a way that the sub-vector of β̂∗(α) indexed by α is equal to β̂(α). Then it is easy to see that proving (7) is equivalent to proving lim inf n→∞ n−1H(β̂∗(α), n)> 0 a.s. for any incorrect model α ∈Aw . (54) Define A0 = {β : ||β − β0|| ≤ 1 2 min 1≤i≤pα0 |β0(α0)i|= b}. G. Qian / Eur. J. Pure Appl. Math, 3 (2010), 417-434 431 Clearly A0 is a compact set; and for any incorrect model α ∈ Aw we have β̂∗(α) 6∈ A0 because ||β̂∗(α)−β0|| ≥ 2b. Moreover, by Theorem 1, the MLE β̂ for the full model is an interior point of A0 almost surely when n is sufficiently large. Since H(β , n) is convex with respect to β , it follows that H(β̂∗(α), n)≥ inf β∈∂ A0 H(β , n) where ∂A0 is the boundary of A0. Now it is sufficient to prove lim inf n→∞ inf β∈∂ A0 n−1H(β , n)> 0 a.s. (55) in order to prove (54). Using result (iii) of Lemma 1, (8), condition (C.7) and Cauchy-Schwarz inequality one can show that R1(β , n)I(β∈∂A0)≥ 1 2 n ∑ k=1 e−2||xk ||·||β−β0||µ0k[x t k (β − β0)] 2I(β∈∂A0) = 1 2 (β − β0) tX t nMnXn(β − β0)I(β∈∂A0) ≥ 1 2 λ1{X t nMnXn}||β − β0||2I(β∈∂A0)≥ 1 2 b3 b2n, which suggests inf β∈∂ A0 R1(β , n)≥ 1 2 b3 b2n. (56) Following the same line as proving (36) one can show that sup β∈∂ A0 |R2(β , n)| = O( p n log log n) a.s.. (57) Since infβ∈∂ A0 H(β , n) ≥ infβ∈∂ A0 R1(β , n)− supβ∈∂ A0 |R2(β , n)| by (8), it follows from (56) and (57) that (55) is true and consequently (54) is true. 5. Discussion The asymptotic results obtained in this paper are based on the assumption that the re- sponse variable follows a Poisson distribution and the candidate models under consideration for selection are based on those available explanatory variables. In practice, there may exist some lurking variables which also affect the response variable. In this situation, a mixture Poisson distribution may be introduced for modeling the response variable, in which an over- dispersion parameter is used to account for the effects of the lurking variables. We refer to [8, section 6.2.3] for details of the over-dispersion log-linear models. Then the model selec- tion procedure can still focus on those available explanatory variables. Provided that a result similar to that in Lemma 2 can be obtained for the mixture Poisson probability (which may be shown by imitating the proof in [2]), it seems the same results as of Theorems 1 to 4 can also REFERENCES 432 be established for the over-dispersion log-linear models using the same methods employed in this paper. It is also possible to extend our asymptotic results to the Poisson regression models with the link functions other than the log-link, which is still under our investigation but shall be presented somewhere else. In addition to the model selection criteria studied here, there are many other approaches for model selection, e.g. the hierarchical Bayesian approach [5] and LASSO [19] etc., for which our results may not be applicable. Finally, other than determining a model selection criterion and assessing its asymptotic performance, there is a computational issue on how to execute a model selection procedure from many possible candidate models. This becomes especially important when the number of candidate models, often of the magnitude 2p, is enormous. However, a thorough investigation of this issue is beyond the scope of this paper. We refer to [10] and [13] for some results in this area. References [1] H Akaike. Information theory and an extension of the maximum likelihood principle. In B.N. Petrov and F. Csáki, editors, Proceedings of the Second International Symposium on Information Theory, pages 267–281, Budapest, 1973. Akadémia Kiadó. [2] H Bohman. Two inequalities for poisson distributions. Skandinavisk Aktuarietidskrift, 46:47–52, 1963. [3] Y Chow and H Teicher. Probability Theory: independence, interchangeability, martingales. Springer, New York, 3 edition, 1997. [4] E George. Statistics in the 21st Century, chapter The variable selection problem, pages 350–358. Chaptman & Hall/CRC, 2002. [5] E George and R McCullock. Approaches for bayesian variable selection. Statistica Sinica, 7:339–373, 1997. [6] N Johnson and S Kotz. Discrete Distributions. Houghton Mifflin Company, Boston, Mas- sachusetts, 1969. [7] C Mallows. Some comments on cp. Technometrics, 15:661–675, 1973. [8] P McCullagh and J Nelder. Generalized Linear Models. Chapman & Hall, London, 2 edition, 1989. [9] V Petrov. Limit Theorems of Probability Theory: sequences of independent random vari- ables. Oxford University Press, 1995. [10] G Qian. Computations and analysis in robust regression model selection using stochastic complexity. Computational Statistics, 14:293–314, 1999. REFERENCES 433 [11] G Qian and H Künsch. Some notes on rissanen’s stochastic complexity. IEEE Transactions on Information Theory, 44:782–786, 1998. [12] G Qian and Y Wu. Strong limit theorems on model selection in generalized linear re- gression with binomial responses. Statisticsa Sinica, 16:1335–1365, 2006. [13] G Qian and X Zhao. On time series model selection involving many candidate arma models. Computational Statistics & Data Analysis, 51:6180–6196, 2007. [14] C Rao and Y Wo. Model Selection, volume 38 of IMS Lecture Notes - Monograph Series, chapter On model selection (with discussion), pages 1–64. Institute of Mathematical Statistics, Beachwood, Ohio, 2001. [15] C Rao and L Zhao. Linear representation of m-estimates in linear models. The Canadian Journal of Statistics, 20:359–368, 1992. [16] J Rissanen. Stochastic Complexity in Statistical Inquiry. World Scientific Publishing Co. Pte. Ltd., Singapore, 1989. [17] J Rissanen. Fisher information and stochastic complexity. IEEE Transactions Information Theory, 42:40–47, 1996. [18] G Schwarz. Estimating the dimension of a model. Annals of Statistics, 6:461–464, 1978. [19] R Tibshirani. Regression shrinkage and selection via the LASSO. Journal of the Royal Statistical Society. [20] Y Wu and M Zen. A strongly consistent linear model selection procedure based on m-estimation. Probability Theory and Related Fields, 113:599–625, 1999. Appendix Proof. (Lemma 1) Since K(t, s) = et − es − es(t − s), it is easy to see that K(s, s) = 0, K ′t(t, s) = et − es, K ′t(s, s) = 0 and K ′′t (t, s) = et > 0. Hence, K(t, s) is strictly convex with respect to t, and K(t, s) ≥ 0 with K(t, s) = 0 only if t = s. Now let F(t, s) = K(t, s)− 1 2 es−2∆(t − s)2 G(t, s) = K(t, s)− 1 2 es+2∆(t − s)2 REFERENCES 434 For any real numbers s, t and ∆ > 0, suppose |t − s| ≤∆. Then we have s− 2∆ < s−∆ ≤ t ≤ s+∆ < s+ 2∆. It is clear that F(s, s) = 0, F ′t(t, s) = et − es − es−2∆(t − s), F ′t(s, s) = 0 and F ′′t (t, s) = et − es−2∆ > 0 Therefore, F(t, s) ≥ 0 with F(t, s) = 0 only if t = s, namely K(t, s) ≥ 1 2 es−2∆(t − s)2. Similarly, we can show that G(t, s) ≤ 0 with G(t, s) = 0 only if t = s, namely K(t, s) ≤ 1 2 es+2∆(t − s)2. This concludes the proof.