Model comparison for generalized linear models with dependent observations

0.31 MB·

by  Shoichi Eguchi

MODEL COMPARISON FOR DEPENDENT GENERALIZED LINEAR MODEL
SHOICHIEGUCHI
Abstract. The classical Bayesian information criterion (BIC) is derived through the stochastic ex-
pansion of marginal likelihood function under suitable regularity condition when models are correctly
specified. However, despite of its popularity, mathematical validity of BIC for possibly misspecified
models with complicated dependence structure is often ignored. Thus it is important to extend the
6 reachoftheclassicalBICwithrigoroustheoreticalfoundationwithallowingmodelmisspecificationand
1 asymptotic normality of estimator. In this paper, we will prove the stochastic expansion of marginal
0 quasi-likelihood function associated with a class of possibly misspecified generalized linear models for
2 dependent data.

v
o 1. Introduction
N
Generalized linear model (GLM, McCullagh and Nelder [25]) is an extension of a linear regression
7
1 model. This model depends on canonical parameter and dispersion parameter, where the former is
representedbythelinkfunctiondeterminedbytheconditionaldistributionoftheresponsevariablegiven

]
T the explanatory variable. Moreover, GLM has many applications and extensions; for example, actuarial
S science(AntonioandBeirlant[3],HabermanandRenshaw[17]),GLMixedMinriskmanagement(McNeil
.
h and Wendin [26]) and generalized additive models (Berg [5], Hastie and Tibshirani [18]).
at We consider data (yj,xj)nj=1 = (yj,xj,1,...,xj,p)nj=1, where yj’s and xj’s are realizations of the re-
m sponse variables Y = (Y ,...,Y ) and the explanatory variables X = (X ,...,X ), respectively,
n 1 n ′ n 1 n ′
[ where the notation means the transpose. Furthermore, we will assume that the conditional distribu-
′

2 tion of Yn given Xn is given by a GLM. Then the conditional distribution is assumed to belong to an
v
exponential family, for example normal, binomial, Poissonand so on. In this paper, we will give a result
2
8 about the stochastic expansion the stochastic expansion of marginalquasi-likelihoodfunction associated
0 with a class of possibly misspecified GLMs for dependent data. Based on the expansion, we propose the
1
quasi-Bayesianinformationcriterion, which is an extension of the generalizedBIC given by Luv and Liu
0
. [23].
1
0 Suppose that we are given M Bayesian candidate models M1,...,MM. Each Mm is described by
6 p ,π (θ),H (θ) θ Θ , where p is the non-zero prior relative occurrence probability of mth-
m m m,n m m
1 ∈
: (cid:8)m(cid:0)odelamongtheM B(cid:1)a(cid:12)yesianm(cid:9)odels,πm istheprior-probabilitydensityonΘm andHm,nisthelogarith-
v (cid:12)
mic quasi-likelihood function. Here “quasi” means that we consider the parametric model that may not
i
X
include the true distribution ofobservationdata. The conventionalBayesianprinciple of model selection
r
a for M1,...,MM is to choose the model that is most likely in terms of the posterior probability, i.e. to
choose model M such that m =argmax P(M y ), where
m0 0 m∈{1,...,M} m| n
exp H (θ) π (θ)dθ p
P(M y )= Θm { m,n } m m ,
m| n (cid:16)MR exp H (θ) π (θ)d(cid:17)θ p
i=1 Θi { i,n } i i
where exp H (θ) π (θ)dθ is calledPthe m(cid:16)aRrginal quasi-likelihood f(cid:17)unction. When the prior plau-
Θm { m,n } m
sibilities on the M competing models would be equal, we select the model that maximizes the marginal
R
quasi-likelihoodfunction;evenifthepriorprobabilitiesarenotequal,wecantriviallycorrecttheselection
manner by the factors p . Hence we focus on the logarithm of the marginal quasi-likelihood function
m
log exp H (θ) π (θ)dθ
m,n m
{ }
(cid:18)ZΘ (cid:19)
Date:November18,2016.
Keywords andphrases. AsymptoticBayesianmodelcomparison,quasi-likelihood,dependentdata,modelmisspecifica-
tion,generalizedlinearmodel.
1
2 SHOICHIEGUCHI
to make model selection.
As was explained in [23], another interpretation of model selection is possible through the Kullback-
Leibler divergence(KL divergence). The KL divergencebetween the true conditionnalmodel g and the
n
marginal quasi-likelihood function exp H (θ) π (θ)dθ is given by
Θ { m,n } m
R
I g ; exp H (θ) π (θ)dθ) =E[logg (Y X )]+E exp H (θ) π (θ)dθ , (1)
n m,n m n n n m,n m
{ } | − { }
(cid:18) ZΘ (cid:19) (cid:20) ZΘ (cid:21)
where the expectation is taken with respect to the true distribution G . Because of (1), we see that
n
exp H (θ) π (θ)dθ is an unbiased estimator of I g ; exp H (θ; ) π(θ)dθ except for a con-
− Θ { m,n } m n Θ { n · }
stant term free of θ. Note that (1) holds true regardless of whether or not the true model is in the set
R (cid:0) R (cid:1)
of candidate models, implying that Bayesian principle of model selection can be restated as choosing
the model that minimizes the KL divergence of the marginal quasi-likelihood function from the true
distribution.
In particular, assume that X is absent and that H (θ) = n logf (y ;θ) for the case of
n m,n j=1 m,n j
independent observations y ,...,y with correctly specified regular models, then Schwarz [27] showed
1 n P
that the marginal quasi-likelihood, log( exp H (θ) π (θ)dθ) admits the stochastic expansion
Θ { m,n } m
R n
p
log exp H (θ) π (θ)dθ = logf (y ;θˆMLE) logn+O (1), (2)
{ m,n } m m,n j m,n − 2 p
(cid:18)ZΘ (cid:19) j=1
X
with θˆMLE denoting the maximum likelihood estimator of θ, under some regularity conditions. Due to
m,n
(2), we obtain the classical Bayesian information criterion for model selection:
n
BIC= 2 logf (y ;θˆMLE)+plogn.
− m,n j m,n
j=1
X
In the past, many authors have investigated the information criteria for model selection in various
settings;see,forexample,BurnhamandAnderson[7]foranaccountofthesedevelopments. Bozdogan[6]
showedthatAkaikeinformationcriterion(AIC,Akaike[1],[2])hasapositiveprobabilityofoverestimating
the true dimension. Casella et al. [8] and Fasen and Kimmig [16] as well as the references therein
studied the model selection consistency of BIC. Moreover, various extensions of AIC and BIC have
been introduced; for example, the extended BIC for large model spaces (Chen and Chen [10]), the
generalized information criteria (Konishi and Kitagawa [21]), the generalized BIC in misspecified GLMs
for independent data (Lv and Liu [23]) and the information criteria in the case of dependent data (e.g.
Sei and Komaki [28] and Uchida [29]).
The rest of the paper is organized as follows. In Section 2, we describe our working model, notations
and assumptions. We also discuss the asymptotic properties of the quasi-maximum likelihood estimator
in possibly misspecified GLM for dependent data. Section 3 presents the stochastic expansion of the
logarithmic marginal quasi-likelihood and the consistency of the model selection with respect to the
optimal model. In Section 4, we illustrate the performance of model selection criterion in correctly
specified and misspecified models. Section 5 presents a real data example. The proofs of our results are
given in Section 6.
2. Quasi-maximum likelihood estimation of dependent GLM
Let Y = (Y ,...,Y ) be the n-dimensional random vector and X = (X ,...,X ) be the n p
n 1 n ′ n 1 n ′
×
random time series. We write X = (X ,...,X ) for any j. We assume that the unknown true
j j,1 j,p ′
distribution of (X ,Y ) has the density g with respect to some dominating σ-finite measure:
n n n
g (x ,y )=g (x )g (y x ),
n n n n n n n n
|
where x =(x ,...,x ), x =(x ,...,x ) and y =(y ,...,y ).
n 1 n ′ j j,1 j,p ′ n 1 n ′
MODEL COMPARISON FOR DEPENDENT GENERALIZED LINEAR MODEL 3
2.1. Model setup. We consider possibly misspecified M candidate models to estimate the true model
G . Each candidate model is given by
n
n
f (x ,y ;θ)=f (x )f (y x ;θ)=f (x ) f (y x ;θ) (3)
m,n n n n n m,n n n n n m,n,j j j
| |
j=1
Y
with θ = (θ1,...,θpm) ∈ Θm, where the mth parameter space Θm ⊂ Rpm is a bounded convex domain
and p p. Although the true structure is unknown, (3) means that the candidate models are given as
m
≤
if Y ,...,Y are (X ,...,X )-conditionally independent and each (X ,...,X )-conditional distribution
1 n 1 n 1 n
ofY depends ononlyX . InapplicationsofGLM,the settingthatY ’sareindependentorconditionally
j j j
independentisoftenassumed. Thelaterconditionisusedtosimplify theoreticalconsiderationandisnot
essential. Furthermore, since X is a multidimensional variable, Y can be given by various expressions.
j j
(3)alsomeansthatwedonotfocusonthedistributionofX ,andthesimilarsituationmaybeconsidered
n
in applications. Therefore, we consider only the true conditional distribution of Y given X and use
n n
GLM M as our working model, with respect to some dominating measure: M is given by
m m
n n
f (y x ;θ)= f (y x ;θ)= exp y x θ b (x θ)+c (y ) , (4)
m,n n| n m,n,j j| j j ′j − m ′j m j
j=1 j=1
Y Y (cid:0) (cid:1)
where, for brevity, we write x θ = pm x θ with d (m),...,d (m) 1,...,p for any m,
′j i=1 j,di(m) i { 1 pm } ⊂ { }
b () and c () are determined by each assumed conditional distribution of Y given X and b () is
m · m · P n n m ·
a sufficiently smooth convex function defined on R; for example, b (θ) = θ2/2 (Gaussian regression)
m
and b (θ) = log(1+eθ) (Logistic regression). For simplicity, we assume that b (θ) = = b (θ) and
m 1 M
···
c (y) = = c (y), that is, only consider variable selection concerning X . As is well known, for any
1 M n
···
n-dimensional random vector Z whose conditional distribution given X is (4), we have E [Z X ] =
n n θ j j
|
∂b (X θ) and V [Z X ]=∂2b (X θ), where ∂b (x)= ∂ b (θ) .
m j′ θ j| j m j′ m ∂θ m θ=x
Since any candidate model M is possibly misspecified and c () of (4) is independent of θ, we may
m m(cid:12) ·
and do define the logarithmic marginal quasi-likelihood function(cid:12)H by
m,n
n
H (θ)= Y X θ b (X θ) . (5)
m,n j j′ − m j′
j=1
X(cid:0) (cid:1)
Any random mapping θˆ such that
m,n
θˆ argmaxH (θ)
m,n m,n
∈
θ∈Θm
is called the quasi-maximum likelihood estimator (QMLE) associated with H . Clearly, when b is
m,n
differentiable, θˆ is the solution to the quasi-score function
m,n
n
∂ H (θ)= Y ∂b (X θ) X =0,
θ m,n j − m j′ j
j=1
X(cid:0) (cid:1)
where ∂ =∂/∂θ.
θ
For notational brevity, from now on we will omit the model index “m” from the notation.
2.2. AsymptoticbehavioroftheQMLE. Inthissection,wewillshowthattheasymptoticproperties
of the QMLE in misspecified GLMs with the dependent observations. Fahrmeir and Kaufmann [15]
studied the consistency and asymptotic normality of the MLE in correctly specified GLMs. Moreover,
White[31]andDomowitzandWhite[12]investigatedthepropertiesoftheQMLEinmisspecifiedmodels
and treated independent data and dependent data, respectively. The settings of Domowitz and White
[12]aremoregeneralthanoursettings,butwewillassumeclearerconditionsmoretailoredfortheGLMs.
In this paper, we will focus on not only the QMLE but also the model selection.
Let =σ(Y ,X ;i j) denote the σ-field representing the data information at stage j. If a and b
j i i n n
F ≤
satisfy a Cb for some constant C >0, we write a .b . We introduce the following conditions:
n n n n
≤
Assumption 2.1. For some constant C 0 and C 0, (i) max ∂ib(x) .1+ xC,
′
≥ ≥ i 1,2,3 | | | |
(ii) E Y 3 σ(X ) .1+ X C′ a.s. for any j N,∈{ }
j j 1 j j
| | |F − ∨ | | ∈
(cid:2) (cid:3)
4 SHOICHIEGUCHI
(iii) supE X 3C+C′+3 < .
j
j N | | ∞
∈
(cid:2) (cid:3)
Assumption2.2. ThereexistsameasurablefunctionF :Rp RsuchthatE[Y σ(X )]=F(X )
j j 1 j j
→ |F − ∨
for every j N.
∈
Assumption 2.3. Denote ζ =(X ,Y ) for any j. For some c>0,
j j j
α(k) c 1e ck
− −
≤
for all k N, where
∈
α(k):=sup sup P[A B] P[A]P[B].
j∈NBA∈∈σ(σζ(iζ;ii≥;i≤j+j)k)| ∩ − |
When Assumption 2.3 holds, ζ ;j =1,2,... is called exponential α-mixing. In particular, Assump-
j
{ }
tion 2.3 implies that ψ := Y F(X ) X , j N, is exponential α-mixing.
j j j j
− ∈
Assumption 2.4. There e(cid:0)xists a non-d(cid:1)egenerate probability measure ν such that the following holds:
n
1
P
(i) F(X )X θ b(X θ) F(x)xθ b(xθ) ν(dx) 0 for any θ Θ,
n j j′ − j′ − ′ − ′ −→ ∈
j=1 Z
X(cid:0) (cid:1) (cid:0) (cid:1)
n
1
(ii) n ∂2b(Xj′θ)XjXj′ − ∂2b(x′θ)xx′ν(dx)−P→0 for any θ ∈Θ,
j=1 Z
X
P
where the notation means the convergence in probability.
−→
Assumption 2.5. There exists a function b:Rp (0, ), (i) for any x, inf ∂2b(xθ) b(x),
′
→ ∞ θ Θ ≥
n ∈
1
(ii) for some constant λ > 0, limsupP λ b(X )X X < λ = 0, where λ () denotes
0 n min n j j j′ 0 min ·
→∞ (cid:20) (cid:18) Xj=1 (cid:19) (cid:21)
the smallest eigenvalues of a given matrix.
Assumptions 2.1, 2.2 and 2.3 are for establishing the asymptotic properties such as Lemmas 6.1 and
6.2. Moreovertheseassumptionsfacilitatethederivationoftheconsistencyandtheasymptoticnormality
of the QMLE θˆ , which are given in Theorems 2.8 and 2.10, respectively.
n
If Assumptions 2.1 and 2.2 hold, for some constant C 0 and C 0,
′
≥ ≥
n n
1 1
supE sup ∂ F(X )X θ ∂b(X θ) sup E F(X ) +sup ∂b(X θ) X
 (cid:12) θn j j′ − j′ (cid:12)≤ n j j′ | j|
n>0 θ∈Θ(cid:12)(cid:12)  Xj=1(cid:0) (cid:1)(cid:12)(cid:12) n>0 Xj=1 (cid:20)(cid:16)(cid:12) (cid:12) θ∈Θ(cid:12) (cid:12)(cid:17) (cid:21)
 (cid:12) (cid:12) n (cid:12) (cid:12) (cid:12) (cid:12)
(cid:12)(cid:12)  (cid:12)(cid:12) .sup 1 E (1+ Xj )C′ +(1+ Xj )C Xj
n | | | | | |
n>0 Xj=1 (cid:20)(cid:16) (cid:17) (cid:21)
< . (6)
∞
Since (6) gives the tightness of 1 n F(X )X θ ∂b(X θ) , we have
n j=1 j j′ − j′
n (cid:8) P (cid:0) (cid:1)(cid:9)
1
P
sup F(X )X θ b(X θ) F(x)xθ b(xθ) ν(dx) 0 (7)
(cid:12)n j j′ − j′ − ′ − ′ (cid:12)−→
θ∈Θ(cid:12)(cid:12) Xj=1(cid:0) (cid:1) Z (cid:0) (cid:1) (cid:12)(cid:12)
(cid:12) (cid:12)
under Assumptions 2.(cid:12)1, 2.2 and 2.4 (i). By using (7), we can prove Theorem 2.8(cid:12)easily.
(cid:12) (cid:12)
From Assumptions 2.1-2.3 and 2.4 (i), we have
n n
1 1 1
nHn(θ)= n ψj′θ+ n F(Xj)Xj′θ−b(Xj′θ)
j=1 j=1
X X(cid:8)(cid:0) (cid:1)(cid:9)
n
1 1
=O + F(X )X θ b(X θ) P F(x)xθ b(xθ) ν(dx)=:H (θ). (8)
p √n n j j′ − j′ −→ ′ − ′ 0
(cid:18) (cid:19) j=1 Z
X(cid:8)(cid:0) (cid:1)(cid:9) (cid:0) (cid:1)
The proof of the tightness of 1 n ψ is given in Lemma 6.1.
√n j=1 j
(cid:8) P (cid:9)
MODEL COMPARISON FOR DEPENDENT GENERALIZED LINEAR MODEL 5
Assumption 2.5 ensures the uniqueness of the QMLE θˆ by the strict concavity of H . Since As-
n n
sumption 2.4 (ii) gives 1∂2H (θ) = ∂2H (θ)+o (1) for any θ, Assumptions 2.1-2.5 imply that the
−n θ n − θ 0 p
equation
∂ H (θ)= F(x) ∂b(xθ) xν(dx)=0
θ 0 ′
−
Z
(cid:0) (cid:1)
admits a unique solution. Then we may define the optimal parameter θ as the unique maximizer of
0
H (θ):
0
θ =argmaxH (θ)
0 0
{ }
θ Θ
∈
The quasi-observed information is given by Γ := 1∂2H (θ ) = 1 n ∂2b(X θ )X X , so that Γ
n −n θ n 0 n j=1 j′ 0 j j′ n
satisfies the equation
P
Γ =Γ +o (1),
n 0 p
where Γ := ∂2b(xθ )xxν(dx).
0 ′ 0 ′
R
Remark 2.6. The β-mixing coefficients of ζ are defined by
j
{ }
β(k):=supE sup P B σ(ζ ;i j) P(B) .
i
j∈N (cid:20)B∈σ(ζi;i≥j+k)(cid:12) (cid:0) | ≤ (cid:1)− (cid:12)(cid:21)
If β(k) = O(e ak) for some a > 0 and for all k (cid:12) N, then ζ is called e(cid:12)xponential β-mixing (e.g.
− j
∈ { }
Davydov [11] and Liebscher [22]). The exponential β-mixing property implies the exponential α-mixing
property. When we replace Assumption 2.3 by the condition that ζ is exponential β-mixing under
j
{ }
someappropriatemomentcondition,thefollowingconditionsfollowonapplyinganobviousdiscrete-time
counterpart of Masuda [24, Lemma 4.3]: (i) For some constant β >0 and q >0,
1 1
1 n q1
supE nβ1sup F(X )X θ b(X θ) F(x)xθ b(xθ) ν(dx) < .
n>0 (cid:20)(cid:18) θ∈Θ(cid:12)(cid:12)nXj=1(cid:0) j j′ − j′ (cid:1)−Z (cid:0) ′ − ′ (cid:1) (cid:12)(cid:12)(cid:19) (cid:21) ∞
(cid:12) (cid:12)
(ii) For some constant β2(cid:12)>0 and q2 >0, (cid:12)
1 n q2
supE nβ2sup ∂2b(X θ)X X ∂2b(xθ)xxν(dx) < .
n j′ j j′ − ′ ′ ∞
n>0 (cid:20)(cid:18) θ∈Θ(cid:12)(cid:12) Xj=1 Z (cid:12)(cid:12)(cid:19) (cid:21)
(cid:12) (cid:12)
Because of Borel-Canntellilemm(cid:12)a, if q1 and q2 can be taken large enough, w(cid:12)e may deduce almost surely
n
1
sup F(X )X θ b(X θ) F(x)xθ b(xθ) ν(dx) 0
n j j′ − j′ − ′ − ′ →
θ∈Θ(cid:12)(cid:12) Xj=1(cid:0) (cid:1) Z (cid:0) (cid:1) (cid:12)(cid:12)
(cid:12) (cid:12)
and (cid:12) (cid:12)
n
1
sup ∂2b(X θ)X X ∂2b(xθ)xxν(dx) 0.
n j′ j j′ − ′ ′ →
θ∈Θ(cid:12)(cid:12) Xj=1 Z (cid:12)(cid:12)
(cid:12) (cid:12)
(cid:12) (cid:12) (cid:3)
Remark 2.7. WecanrelaxAssumption2.2byreplacingE[Y σ(X )]=F(X )withE[Y
j j 1 j j j j 1
|F − ∨ |F − ∨
σ(X )] = F(X ,...,X ) for some m 1. When we assume relaxed condition, Assumption 2.4 (i)
j j m+1 j
− ≥
is modified. For example, in the case of m=2, Assumption 2.4 (i) is given by
n
1
P
F(X ,X )X θ b(X θ) F(x ,x )x θ b(x θ) P(x ,dx )ν(dx ) 0,
n j−1 j j′ − j′ − 1 2 ′2 − ′2 1 2 1 −→
j=2 Z
X(cid:0) (cid:1) (cid:0) (cid:1)
where P is a transition function. (cid:3)
Theorem 2.8. Under Assumptions 2.1-2.5, the QMLE satisfies
θˆ P θ
n 0
−→
as n .
→∞
6 SHOICHIEGUCHI
Assumption 2.9. (i) X ;j =1,2,... is strictly stationary.
j
{ }
n n
1 ′
(ii) For some Σ >0, E Y ∂b(X θ ) X Y ∂b(X θ ) X Σ .
0 n j − j′ 0 j j − j′ 0 j → 0
(cid:20)(cid:26)j=1 (cid:27)(cid:26)j=1 (cid:27)(cid:21)
X(cid:0) (cid:1) X(cid:0) (cid:1)
Theorem2.10. UnderAssumptions2.1-2.5and2.9,theasymptoticdistributionoftheQMLEisnormal:
√n(θˆn−θ0)−→L N(0,Γ−01Σ0Γ−01).
When the candidate model is correctlyspecified, Σ0 =Γ0, i.e. √n(θˆn−θ0)−→L N(0,Γ−01). The reason
is because the correctly specified model gives the equation
n n
1 ′
E Y ∂b(X θ ) X Y ∂b(X θ ) X
n j − j′ 0 j j − j′ 0 j
(cid:20)(cid:26)j=1 (cid:27)(cid:26)j=1 (cid:27)(cid:21)
X(cid:0) (cid:1) X(cid:0) (cid:1)
n
1 2
= E E Y ∂b(X θ ) σ(X ) X X
n j − j′ 0 |Fj−1∨ j j j′
Xj=1 h (cid:2)(cid:0) (cid:1) (cid:3) i
1
+2 E Y ∂b(X θ ) E Y ∂b(X θ ) σ(X ) X X
n i− i′ 0 j − j′ 0 |Fj−1∨ j i j′
Xi<j h(cid:0) (cid:1) (cid:2)(cid:0) (cid:1) (cid:3) i
n
1 2
= E E Y E [Y X ] X X X
n θ0 j − θ0 j| j | j j j′
Xj=1 h (cid:2)(cid:0) (cid:1) (cid:3) i
1
+2 nE Yi−∂b(Xi′θ0) Eθ0 Yj −Eθ0[Yj|Xj] |Xj XiXj′
Xi<j h(cid:0) (cid:1) (cid:2)(cid:0) (cid:1) (cid:3) i
n
1
= n E ∂2b(Xj′θ0)XjXj′ +0
Xj=1 h i
= ∂2b(xθ )xxν(dx)
′ 0 ′
Z
under Assumptions 2.1-2.5 and 2.9 (i).
Remark 2.11. Let Assumptions 2.2-2.5 and 2.9 (i) hold. The condition
n
1 2
E Y ∂b(X θ ) X X Σ 0
(cid:12)n j − j′ 0 j j′ − 0(cid:12)→
(cid:12) j=1 (cid:12)
(cid:12) X (cid:2)(cid:0) (cid:1) (cid:3) (cid:12)
(cid:12) (cid:12)
implies Assumption 2.9 (ii), si(cid:12)nce it follows from this condition and D(cid:12)oukhan [13, Theorem 3] that
(cid:12) (cid:12)
n n
1 ′
E Y ∂b(X θ ) X Y ∂b(X θ ) X Σ
(cid:12)n j − j′ 0 j j − j′ 0 j − 0(cid:12)
(cid:12) (cid:20)(cid:26)j=1 (cid:27)(cid:26)j=1 (cid:27)(cid:21) (cid:12)
(cid:12) X(cid:0) (cid:1) X(cid:0) (cid:1) (cid:12)
(cid:12) (cid:12)
(cid:12)(cid:12) 1 n 2 2 (cid:12)(cid:12)
E Y ∂b(X θ ) X X Σ + cov Y ∂b(X θ ) X ,(Y ∂b(X θ ) X
≤(cid:12)n j − j′ 0 j j′ − 0(cid:12) n i− i′ 0 i j − j′ 0 j
(cid:12)(cid:12)(cid:12) Xj=1 (cid:2)(cid:0) (cid:1) (cid:3) (cid:12)(cid:12)(cid:12) Xi<j(cid:12)(cid:12) (cid:2)(cid:0) (cid:1) (cid:1) (cid:3)(cid:12)(cid:12)
(cid:12)(cid:12)1 n 2 (cid:12)(cid:12) 2 (cid:12) (cid:12)
. E Y ∂b(X θ ) X X Σ + α(j i) 0
(cid:12)n j − j′ 0 j j′ − 0(cid:12) n − →
(cid:12) j=1 (cid:12) i<j
(cid:12) X (cid:2)(cid:0) (cid:1) (cid:3) (cid:12) X
as n (cid:12)(cid:12) . (cid:12)(cid:12) (cid:3)
→(cid:12)∞ (cid:12)
3. Quasi-Bayesian information criterion for dependent GLM
3.1. Stochastic expansion. We use the GLM as our working model to choose the optimal model, so
we consider the stochastic expansion of the marginal quasi-likelihood in GLM.
n
1
Assumption 3.1. F(X ) ∂b(X θ ) X =O (1).
√n j − j′ 0 j p
j=1
X(cid:0) (cid:1)
The next theorem shows the asymptotic behavior of the log marginal quasi-likelihood function.
MODEL COMPARISON FOR DEPENDENT GENERALIZED LINEAR MODEL 7
Theorem 3.2. Assume that Assumptions 2.1-2.5 and 3.1 hold and that the following conditions are
satisfied:
(i) π(θ )>0, supπ(θ)< .
0
∞
θ Θ
∈
u
(ii) For every M >0, sup π θ + π(θ ) 0 as n .
0 0
√n − → →∞
(iii) logπ(θˆn) logπ(θ0|u)|<=Mop(cid:12)(cid:12)(cid:12)(1(cid:18)). (cid:19) (cid:12)(cid:12)(cid:12)
− (cid:12) (cid:12)
Then we have the expansion
log exp H (θ) π(θ)dθ
n
{ }
(cid:18)ZΘ (cid:19)
n n
p p 1 1
= Y X θˆ b(X θˆ ) logn+ log2π logdet ∂2b(X θˆ )X X +logπ(θˆ )+o (1)
j j′ n− j′ n − 2 2 − 2 n j′ n j j′ n p
j=1 (cid:18) j=1 (cid:19)
X(cid:0) (cid:1) X
n n
p 1
= Y X θˆ b(X θˆ ) + log2π logdet ∂2b(X θˆ )X X +logπ(θˆ )+o (1).
j j′ n− j′ n 2 − 2 j′ n j j′ n p
j=1 (cid:18)j=1 (cid:19)
X(cid:0) (cid:1) X
Remark 3.3. Suppose that we replace Assumptions 2.4 and 2.5 (ii) by the following conditions:
n
1
(i) F(X )X θ b(X θ) F(x)xθ b(xθ) ν(dx) almost surely as n , uniformly in
n j j′ − j′ → ′ − ′ →∞
j=1 Z
X(cid:0) (cid:1) (cid:0) (cid:1)
θ Θ.
∈
n
1
(ii) ∂2b(X θ)X X ∂2b(xθ)xxν(dx) almost surely as n , uniformly in θ Θ.
n j′ j j′ → ′ ′ →∞ ∈
j=1 Z
X
n
1
(iii) For some constant λ >0, P limsupλ b(X )X X <λ =0.
0 min n j j j′ 0
n
(cid:20) →∞ (cid:18) Xj=1 (cid:19) (cid:21)
Then, as in Cavanaugh and Neath [9, Section 3], we can show that the log marginal quasi-likelihood
function almost surely satisfies the expansion similar to Theorem 3.2, i.e. almost surely
n
p p
log exp{Hn(θ)}π(θ)dθ = YjXj′θˆn−b(Xj′θˆn) − 2logn+ 2log2π
(cid:18)ZΘ (cid:19) j=1
X(cid:0) (cid:1)
n
1 1
logdet ∂2b(X θˆ )X X +logπ(θˆ )+o(1).
− 2 n j′ n j j′ n
(cid:18) j=1 (cid:19)
X
(cid:3)
DuetoTheorem3.2,wedefinethequasi-Bayesianinformationcriterion(QBIC)andBICfordependent
GLM by
n n
QBIC=−2 YjXj′θˆn−b(Xj′θˆn) +logdet ∂2b(Xj′θˆn)XjXj′
j=1 (cid:18)j=1 (cid:19)
X(cid:0) (cid:1) X
and
n
BIC= 2 Y X θˆ b(X θˆ ) +plogn.
− j j′ n− j′ n
j=1
X(cid:0) (cid:1)
Let QBIC(1),...,QBIC(M) be the QBIC for each candidate model. We calculate QBIC(1),...,QBIC(M)
and select the best model M having the minimum-QBIC value:
m0
m = argmin QBIC(m).
0
m 1,...,M
∈{ }
We can also select the best model by using BIC in a similar manner. As directly seen by the definition,
theQBIChavemorecomputationalloadthantheBIC.SincetheQBICinvolvestheobserved-information
matrixquantity,whichisdirectlycomputedfromdata,theQBICwouldmoreeffectivelytakedatadepen-
dence into account. Furthermore,the penalty (second-term)ofQBIC consists ofthe secondderivativeof
8 SHOICHIEGUCHI
H with respectto θ, sothe dimensionofthe parameteraffectsthe penalty. Thatis,the QBICconsiders
n
the complexity of the model.
3.2. Model selection consistency. Let Θi Rpi and Θj Rpj be the parameter space associated
⊂ ⊂
with Mi and Mj, respectively. If pi <pj and there exist a matrix A∈Rpj×pi with A′A=Ipi×pi as well
as a c Rpj such that Hi,n(θ)=Hj,n(Aθ+c) for all θ Θi, we say that Θi is nested in Θj.
∈ ∈
Under Assumptions 2.1-2.5, when m satisfies
0
m =argmindim(Θ ),
0 m
{ }
m
∈M
where =argmax H (θ )=argmax F(x)xθ b (xθ ) ν(dx), wesay
M m∈{1,...,M} m,0 m,0 m∈{1,...,M} ′ m,0− m ′ m,0
that M is the optimal model.
m0 R (cid:0) (cid:1)
Theorem 3.4. Assume that Assumptions 2.1-2.5, 3.1 are satisfied and that there exists a unique m
0
∈
1,...,M such that M is the optimal model. For any fixed m 1,...,M m , if Θ is nested
{ } m0 ∈ { }\{ 0} m0
in Θ , or H (θ)=H (θ ) for any θ Θ , then
m m,0 6 m0,0 m0,0 ∈ m
lim P QBIC(m0) QBIC(m) <0 =1.
n −
→∞
(cid:2) (cid:3)
This theorem implies that the probability that the optimal model is selected by using QBIC tends to
1 as n . The probability that BIC choose the optimal model can be handled analogously.
→∞
4. Examples and Simulation results
In this section, we conduct simulations to evaluate finite sample performance of the model selection
by using QBIC, BIC and formal AIC (fAIC). Since we do not deal with the theoretical part of AIC in
this paper, we use the word fAIC as AIC, i.e. fAIC of mth model is defined by
fAIC(m) = 2H (θˆ )+2p .
m,n m,n m
−
Let θ be the true value. We here set the initial value as the value generated from uniform distribution
∗
U(θ 1, θ +1) to use optimat softwareR for numerical optimization. Note that it may happens that
∗ ∗
−
θ =θ when the model is misspecified.
0 ∗
6
Remark 4.1. Since the true model is unknown in practice, we can not explicitly compute the optimal
parameter and the optimal model. However, the QMLE, QBIC and BIC are computable by assuming
specific candidate models. (cid:3)
4.1. Model selection in correctly specified model. We assume that the explanatory variables
X ,...,X are given by
j,1 j,4
X =1(j 1),
j,1
≥
X =1, X =0.5X +ǫ , (j 2),
1,2 j,2 j 1,2 j,2
− ≥
X =0, X = 0.7X +ǫ , (j 2),
1,3 j,3 j 1,3 j,3
− − ≥
X = 1, X =0.8X +ǫ , (j 2),
1,4 j,4 j 1,4 j,4
− − ≥
where the error vector (ǫ ,ǫ ,ǫ ) N(0,Σ) with Σ = (0.5k ℓ) . Moreover, the response
j,2 j,3 j,4 | − | k,ℓ=1,2,3
∼
variable Y is obtained from the true model defined by the linear logistic regressionmodel
j
exp(X θ )
j′ ∗
Y B 1, , (9)
j
∼ 1+exp(X θ )
(cid:18) j′ ∗ (cid:19)
wherethe true value θ =(0, 3,0,1)andB(1,P)is a Bernoullidistributionwith success probabilityP.
∗
−
We consider this model for the following combination of X :
j
Model1:X =(X ,X ,X ,X ); Model2:X =(X ,X ,X );
j j,1 j,2 j,3 j,4 j j,1 j,2 j,3
Model3:X =(X ,X ,X ); Model4:X =(X ,X ,X ); Model5:X =(X ,X ,X );
j j,1 j,2 j,4 j j,1 j,3 j,4 j j,2 j,3 j,4
Model6:X =(X ,X ); Model7:X =(X ,X ); Model8:X =(X ,X );
j j,1 j,2 j j,1 j,3 j j,1 j,4
MODEL COMPARISON FOR DEPENDENT GENERALIZED LINEAR MODEL 9
Table1. ThenumberofmodelsselectedbyQBIC,BICandfAICinSection4.1over10000simulations
forvariousn(1-15representthemodels,andthetruemodelisModel10)
Criteria n=50
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
∗
QBIC 1489 65 2084 0 1260 201 0 0 56 4666 0 0 172 0 0
BIC 99 18 531 0 562 222 0 0 86 7720 0 0 762 0 0
fAIC 479 55 1310 0 1424 192 0 0 93 6242 0 0 205 0 0
Criteria n=100
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
∗
QBIC 298 2 1483 0 989 10 0 0 4 7206 0 0 8 0 0
BIC 19 1 323 0 397 15 0 0 8 9179 0 0 58 0 0
fAIC 347 1 1380 0 1367 7 0 0 2 6895 0 0 1 0 0
Criteria n=200
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
∗
QBIC 86 0 910 0 616 0 0 0 0 8388 0 0 0 0 0
BIC 5 0 235 0 222 0 0 0 0 9538 0 0 0 0 0
fAIC 281 0 1314 0 1414 0 0 0 0 6991 0 0 0 0 0
Model9:X =(X ,X ); Model10:X =(X ,X ); Model11:X =(X ,X );
j j,2 j,3 j j,2 j,4 j j,3 j,4
Model12:X =X ; Model13:X =X ; Model14:X =X ; Model15:X =X .
j j,1 j j,2 j j,3 j j,4
Then the true model is Model 10, and Models 1, 3, 5 contain the true model. We simulate the number
of the model selected by using QBIC, BIC and fAIC among the candidate Models 1-15 over 10000
simulations. For example, in the case of Model 1, QBIC, BIC and fAIC given by
n 4 4 n exp 4 X θˆ X X
QBIC= 2 Y X θˆ log 1+exp X θˆ +logdet i=1 j,i i j j′
− Xj=1(cid:26) jXi=1 j,i i− (cid:18) (cid:16)Xi=1 j,i i(cid:17)(cid:19)(cid:27) Xj=1 1+(cid:16)exPp 4i=1X(cid:17)j,iθˆi 2
 
n 4 4  (cid:16) (cid:16)P (cid:17)(cid:17) 
BIC= 2 Y X θˆ log 1+exp X θˆ +4logn,
j j,i i j,i i
− −
Xj=1(cid:26) Xi=1 (cid:18) (cid:16)Xi=1 (cid:17)(cid:19)(cid:27)
n 4 4
fAIC= 2 Y X θˆ log 1+exp X θˆ +4 2.
j j,i i j,i i
− − ×
Xj=1(cid:26) Xi=1 (cid:18) (cid:16)Xi=1 (cid:17)(cid:19)(cid:27)
In what follows, we verify the assumptions for QBIC (BIC). In the present situation, the function b
defined in (4) is givenby b(θ)=log(1+eθ), and Assumptions 2.1 and 2.2 are satisfied with E[Y
j j 1
|F − ∨
σM(oXrje)o]v=er,Fth(Xejfu)n=cteio−n3Xbj,o2f+AXjs,s4u/m(1p+tioen−32X.5j,2c+aXnj,b4e).gIivnepnabryticbu(xla)r,=suepjC∈xN/E(1|X+je|Cq x<)2∞forfosroemveercyonqs>tan0t.
− | | (cid:2) (cid:3)| |
C > 0 satisfying sup xθ C x. If (X ,Y );j = 1,2,... has the exponential β-mixing property,
θ∈Θ| ′ | ≤ | | { j j }
Assumptions 2.3, 2.4 and 3.1 hold. Here the exponential β-mixing property of (X ,Y );j =1,2,... is
j j
{ }
shown by Baraud et al. [4] and Doukhan [13, Section 2.4]. The validity of the assumptions of Sections
4.2 and 4.3 can be checked in a similar way.
Table 1 summarizes the comparison results of the frequency of the model selection. Model 10 is
selectedwith high frequency for all criterions andn. Moreover,the probability that Model 10 is selected
by QBIC and BIC becomes higher as n increases. In Table 2, the differences between the true value and
the estimators in specified models are getting small when n gets increased. From these results, we can
observe the consistency of the estimators and the model selection consistency of QBIC and BIC.
4.2. Model selection in misspecified model. We use the same conditions as in the previous section
except for the true model. In this simulation, the response variable Y is obtained from the true model
j
defined by
Yj ∼B 1,Φ(Xj′θ∗) ,
(cid:16) (cid:17)
where Φ(x)= x 1 exp( t2)dt. Then Models 1-15 are misspecified models.
√2π −2
−∞
R
10 SHOICHIEGUCHI
Table 2. Themeanandthestandarddeviation (s.d.) of estimatorsθˆ1,θˆ2,θˆ3 andθˆ4 ineachmodel
forvariousn(1-15representthemodels,andthetrueparameter (θ1∗,θ2∗,θ3∗,θ4∗)=(0,−3,0,1))
n=50 n=100 n=200
θˆ1 θˆ2 θˆ3 θˆ4 θˆ1 θˆ2 θˆ3 θˆ4 θˆ1 θˆ2 θˆ3 θˆ4
1 mean -0.0793 -8.3409 -0.0219 2.7990 0.0004 -3.3727 -0.0057 1.1266 0.0023 -3.1642 0.0004 1.0542
s.d. 8.1943 35.1378 5.9571 12.8200 0.3895 0.8918 0.2889 0.3744 0.2425 0.4946 0.1807 0.2123
2 mean -0.0505 -2.7481 0.1061 – -0.0425 -2.1878 0.0867 – -0.0167 -2.0734 0.0856 –
s.d. 2.0544 7.9339 1.3571 – 0.4653 0.5332 0.2141 – 0.3146 0.3336 0.1395 –
3 mean 0.0355 -5.7372 – 1.8913 0.0001 -3.2941 – 1.0993 0.0021 -3.1332 – 1.0441
s.d. 4.6176 23.3197 – 8.1816 0.3763 0.8210 – 0.3508 0.2397 0.4804 – 0.2078
4 mean -0.0999 – -0.2581 0.3139 -0.0451 – -0.2364 0.2852 -0.0168 – -0.2250 0.2746
s.d. 0.4791 – 0.2318 0.3193 0.3109 – 0.1453 0.1940 0.2100 – 0.0976 0.1264
5 mean – -5.6381 0.0219 1.8650 – -3.2928 -0.0052 1.0989 – -3.1336 0.0003 1.0442
s.d. – 22.5995 3.7461 7.8695 – 0.8333 0.2792 0.3465 – 0.4835 0.1787 0.2057
6 mean -0.0635 -2.4621 – – -0.0429 -2.1293 – – -0.0169 -2.0324 – –
s.d. 1.8364 4.1725 – – 0.4578 0.5074 – – 0.3127 0.3232 – –
7 mean -0.1086 – -0.1960 – -0.0518 – -0.1808 – -0.0217 – -0.1723 –
s.d. 0.4591 – 0.2093 – 0.3188 – 0.1351 – 0.2227 – 0.0921 –
8 mean -0.1006 – – 0.2681 -0.0453 – – 0.2483 -0.0170 – – 0.2415
s.d. 0.4660 – – 0.3022 0.3063 – – 0.1875 0.2081 – – 0.1230
9 mean – -2.3773 0.1058 – – -2.1041 0.0878 – – -2.0342 0.0860 –
s.d. – 5.2112 1.0538 – – 0.5008 0.2061 – – 0.3263 0.1372 –
10∗ mean – -4.2068 – 1.3952 – -3.2211 – 1.0741 – -3.1037 – 1.0344
s.d. – 13.1535 – 4.3787 – 0.7702 – 0.3259 – 0.4699 – 0.2013
11 mean – – -0.2546 0.3124 – – -0.2350 0.2855 – – -0.2243 0.2747
s.d. – – 0.2218 0.2851 – – 0.1424 0.1840 – – 0.0967 0.1230
12 mean -0.0695 – – – -0.0368 – – – -0.0179 – – –
s.d. 0.3186 – – – 0.2548 – – – 0.1936 – – –
13 mean – -2.6475 – – – -2.5688 – – – -2.5278 – –
s.d. – 0.3694 – – – 0.3052 – – – 0.2789 – –
14 mean – – -0.1525 – – – -0.1528 – – – -0.1517 –
s.d. – – 0.1694 – – – 0.1218 – – – 0.0901 –
15 mean – – – 0.5703 – – – 0.5507 – – – 0.5394
s.d. – – – 0.2553 – – – 0.2478 – – – 0.2477
Table3. ThenumberofmodelsselectedbyQBIC,BICandfAICinSection4.2over10000simulations
forvariousn(1-15representthemodels)
Criteria n=100
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
QBIC 965 0 2025 0 1321 0 0 0 0 5689 0 0 0 0 0
BIC 41 0 435 0 398 0 0 0 0 9125 0 0 1 0 0
fAIC 443 0 1452 0 1538 0 0 0 0 6567 0 0 0 0 0
Criteria n=200
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
QBIC 223 0 1338 0 915 0 0 0 0 7524 0 0 0 0 0
BIC 9 0 278 0 274 0 0 0 0 9439 0 0 0 0 0
fAIC 349 0 1436 0 1414 0 0 0 0 6801 0 0 0 0 0
Criteria n=300
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
QBIC 108 0 1009 0 694 0 0 0 0 8189 0 0 0 0 0
BIC 5 0 190 0 216 0 0 0 0 9589 0 0 0 0 0
fAIC 295 0 1352 0 1388 0 0 0 0 6965 0 0 0 0 0
From Table 3, we obtain similar results even though the candidate models do not include the true
model. Table 4 summarizes the mean and the standard deviation of estimators in eachmodel. Since the
optimalparametervalue is not givenhere,we cannot see the differences between the optimal parameter
value and the estimators, although the standard deviations become smaller when n become larger.
4.3. Model selection in univariate time series model. Let X = (Z ,Z ,...,Z ) be the
j j j 1 j (p 1) ′
− − −
explanatory vector for any j 1,...,n , where for every i 2,...,n , Z is given by
i
∈{ } ∈{ }
Z = =Z =0, Z =1, Z =0.6Z +ǫ ,
n+2 0 1 i i 1 i
− ··· −