Bias-corrected AIC for selecting variables in multinomial logistic regression models

13 Pages·2012·0.3 MB·English

by  Hirokazu Yanagihara

LinearAlgebraanditsApplications436(2012)4329–4341
ContentslistsavailableatSciVerseScienceDirect
Linear Algebra and its Applications
journal homepage:www.elsevier.com/locate/laa
Bias-corrected AIC for selecting variables in multinomial
<
logistic regression models
a,∗,1 b a c,2
HirokazuYanagihara ,Ken-ichiKamo ,ShinpeiImori ,KenichiSatoh
aDepartmentofMathematics,GraduateSchoolofScience,HiroshimaUniversity,1-3-1Kagamiyama,Higashi-Hiroshima739-8626,Japan
bDepartmentofLiberalArtsandSciences,SapporoMedicalUniversity,South1,West17,Chuo-ku,Sapporo060-8543,Japan
cDepartmentofEnvironmetricsandBiometrics,ResearchInstituteforRadiationBiologyandMedicine,HiroshimaUniversity,
1-2-3Kasumi,Minami-ku,Hiroshima734-8553,Japan
ARTICLE INFO ABSTRACT
Articlehistory: Inthispaper,weconsiderthebiascorrectionofAkaike’sinformat-
Received22January2012 ion criterion (AIC) for selecting variables in multinomial logistic
Accepted25January2012 regressionmodels.Forsimplifyingaformulaofthebias-corrected
Availableonline24February2012 AIC,wecalculatethebiasoftheAICtoariskfunctionthroughtheex-
SubmittedbyR.A.Brualdi pectationsofpartialderivativesofthenegativelog-likelihoodfunc-
tion. As a result, we can express the bias correction term of the
AMSclassification:
bias-correctedAICwithonlythreematricesconsistingofthesecond,
62H12
third,andfourthderivativesofthenegativelog-likelihoodfunction.
62J12
Byconductingnumericalstudies,weverifythattheproposedbias-
Keywords: correctedAICperformsbetterthanthecrudeAIC.
AIC ©2012ElsevierInc.Allrightsreserved.
Biascorrection
Multinomiallogisticmodel
MLE
Partialdifferentialoperator
Variableselection
1. Introduction
Amultinomiallogisticregressionmodelisaregressionmodelthatgeneralizesalogisticregres-
sion by allowing more than two discrete response variables. When categories are unordered, the
<
Theauthorswouldliketothanktheassociateeditorandthereviewerforvaluablecomments.
∗
Correspondingauthor.
E-mailaddress:[email protected](H.Yanagihara).
1 SupportedbytheMinistryofEducation,Science,Sports,andCulture,Grant-in-AidforChallengingExploratoryResearch,
#22650058,2010–2012.
2 Supported by the Ministry of Education, Science, Sports, and Culture, Grant-in-Aid for Young Scientists (B), #23700337,
2011–2014.
0024-3795/$-seefrontmatter©2012ElsevierInc.Allrightsreserved.
http://dx.doi.org/10.1016/j.laa.2012.01.018
4330 H.Yanagiharaetal./LinearAlgebraanditsApplications436(2012)4329–4341
multinomial logistic model is one strategy often used. The multinomial logistic regression model
hasbeenintroducedinmanytextbooksforappliedstatisticalanalysis(seee.g.[1,Chapter8.1]),and
evennowitiswidelyusedinbiometrics,econometrics,psychometrics,sociometrics,andmanyother
fieldsofapplicationsforthepredictionofprobabilitiesofdifferentpossibleoutcomesofcategorically
distributedresponsevariablesbyasetofexplanatoryvariables(e.g.[2–4]).Inaddition,themultino-
miallogisticregressionmodelcanbeeasilyfittedtoobserveddatabyusingthe“vglm"functionin“R"
[5]thatisafreesoftwareenvironmentforstatisticalcomputingandgraphics.Sincewewouldliketo
specifythefactorsaffectingtheprobabilitiesofresponsevariablesintheregressionanalysis,searching
fortheoptimalsubsetofexplanatoryvariablesisimportant.
Akaike’s information criterion (AIC) proposed by Akaike [6,7] is widely used for selecting the
bestmodelamongthecandidatemodels(fordetailsofstatisticalmodelselection,seee.g.[8–10]).
ThemodelhavingthesmallestAICamongthecandidatemodelsisregardedasthebestmodel.Inthe
multinomiallogisticregressionmodel,thesubsetofexplanatoryvariablesinthebestmodelisthebest
subset.However,theAICmayperformpoorly;thatis,amodelhavingtoomanyparameterstendsto
bechosenasthebestmodelwhenthesamplesizeissmallorthenumberofunknownparametersis
large.Suchaproblemisoftenresolvedbyusingabias-correctedAIC(seee.g.[8,Chapter2.4]).TheAIC
isanestimatoroftheriskfunctionconsistingofpredictiveKullback–Leibler(K–L)information[11],
whichmeasuresthediscrepancybetweenthetruemodelandthecandidatemodel.Theorderofthe
biasoftheAICisO(n−1)whenthecandidatemodelincludesthetruemodel,wherenisthesample
size.AlthoughtheAICisanasymptoticunbiasedestimatoroftheriskfunction,ithasanonnegligible
biastotheriskfunctionwhenthesamplesizeissmallorthenumberofunknownparametersislarge.
Abias-correctedAICcalledCAICinthispaperimprovesthebiasofAICtoO(n−2)undertheassumption
thatthecandidatemodelincludesthetruemodel.
TheCAICinthelogisticregressionmodelswasobtainedbyYanagiharaetal.[12].ButtheCAICin
multinomialregressionmodelshasnotbeenderivedyet,althoughthemultinomiallogisticregression
modeliswidelyusedinmanyapplicationfields.TheCAICcanbeobtainedbyremovingthebiasofthe
AICtotheriskfunctionfromtheAICwiththeuseofaconsistentestimatorofthebias.Thebiasofthe
AICtotheriskfunctionisthenevaluatedbymomentsofthemaximumlikelihoodestimator(MLE)
ofunknownparameters.Sincesuchmomentsarerepresentedbythemomentsofresponsevariables,
calculatingthemomentsofresponsevariablesisessentialforevaluatingthebiasoftheAICintheor-
dinarybiascorrectionmethod,whichisusedin[12,13],etc.However,inthecaseofmultipleresponse
variables,calculationsandexpressionsofthemomentsoftheMLEmediatedbythemomentsofre-
sponsevariablesbecomecomplicated.Hence,withoutdirectlycalculatingthemomentsofresponse
variables,wederivethemomentsoftheMLEbyusingexpectationsofthepartialderivativesofthe
negativelog-likelihoodfunction.Thisdifferentapproachfromtheordinarybiascorrectionmethod
leadstoasimpleexpressionofthebiascorrectiontermoftheCAIC.Infact,thebiascorrectiontermof
ourCAICisrepresentedbyonlythreematricesconsistingofthesecond,third,andfourthderivatives
ofthenegativelog-likelihoodfunction.
Thepresentpaperisorganizedasfollows.InSection2,wegiveastochasticexpansionoftheMLE.
InSection3,theCAICinmultinomiallogisticregressionmodelsisproposed.InSection4,weverify
thattheproposedCAIChasbetterperformancethantheAICbyconductingnumericalexperiments.In
Section5,weconcludeourdiscussion.TechnicaldetailsareprovidedinAppendix.
2. StochasticexpansionofMLE
{ , } ,...,
Supposethatthedataconsistsofasequence y x ,wherey y arer-dimensionalinde-
i i 1 m
,...,
pendent unordered discrete random vectors, and x x are k-dimensional vectors of known
β =(β ,...,β )(cid:2) 1 m
constants.Let beakr-dimensionalunknownregressioncoefficientvectorthatis
β=(β(cid:2)1,...,β(cid:2)k)r(cid:2) β β =(β ,...,
partitionedas 1 r ,where jisak-dimensionalvectordenotedby j (j−1)k+1
β )(cid:2) ( , (cid:2))(cid:2) = ( , ,..., )(cid:2)
dijsktr.i(cid:2)bInuttehdeamccuolrtidnionmgtiaoltlhoegimstiucltriengormesisailodnismtroibduelt,iownewasitshumtheetnhuamt byei0royfieventsyni0i (yni1i = (cid:2)yrj=ir0yiijs,
n= mi=1ni)andthecellprobabilityvector(pi0(β),pi(β)(cid:2))(cid:2),givenby
H.Yanagiharaetal./LinearAlgebraanditsApplications436(2012)4329–4341 4331
1
pi0(β)= 1+(cid:2)rj=1exp(x(cid:2)iβj)
(β)=( (β),..., (β))(cid:2)
pi pi1 pir (1)
(cid:3) (cid:4)
( (cid:2)β ) ( (cid:2)β ) (cid:2)
exp x exp x
= (cid:2) i 1 ,..., (cid:2) i r .
1+ rj=1exp(x(cid:2)iβj) 1+ rj=1exp(x(cid:2)iβj)
β
TheMLEof isobtainedbymaximizingthelog-likelihoodfunction.Inthispaper,weassumethatthe
maximumexists,i.e.,thatthesupremumdoesnotoccurattheboundary.Byomittingtheconstant
term,thelog-likelihoodfunctionofthemultinomiallogisticregressionmodelin(1)isexpressedas
⎧ ⎛ ⎞⎫
(cid:5)m ⎨ (cid:5)r ⎬
(cid:3)(β)= ( ⊗ )(cid:2)β− ⎝ + ( (cid:2)β )⎠ .
y x n log 1 exp x
⎩ i i i i j ⎭
i=1 j=1
β
Hence,theMLEof isgivenby
βˆ = (cid:3)(β).
argmax
β
βˆ
ToevaluateabiasoftheAICtotheriskfunction,astochasticexpansionof isneeded.Thepurpose
ofthissectionistoobtainthestochasticexpansionβˆ uptotheordern−3/2.Twocasesserveasa
frameworkforasymptoticapproximations:
→∞
Case(i):n’sarefixed,andm ,
j
Case(ii):misfixed,n →∞andρ−1 =n/n =O(1)foreachj.
j j j
AlthoughweonlyconsiderCase(i)inthispaper,ourformulacanalsobeappliedtoCase(ii).
Supposethatx ,...,x aremembersofanadmissiblecompactsetF,i.e.,x ,...,x ∈ F.To
1 m 1 m
expandtheMLE,weconsiderthefollowingregularityassumptions(seee.g.[14]):
C.1 :β ∈B,whereBisaconvexandopensetinRk,
C.2 :(I ⊗x)(cid:2)β ∈ (cid:5)0,i = 1,2,...,forallβ ∈ B,where(cid:5)0istheinterioroftheconvexnatural
r i
parameterspace(cid:5)⊂Rr,
∃ . (cid:2) (cid:2) =( ,..., )(cid:2)
C.3 : m st.X Xhasthefullrankform m ,whereX x x .
0 0 1 m
ConditionC.1guaranteestheuniquenessoftheMLEifitexists.ConditionC.2(cid:2)isnecessarytoobtain
pthoesimtivueltdineofimniitaelfloogriasltlicβre∈grBe,smsio(cid:2)nmmod,welhfoerreallβ.ConditionC.3ensuresthat mi=1ni(cid:3)i(β)⊗xix(cid:2)iis
0
(cid:3)(β)=diag(p(β))−p(β)p(β)(cid:2). (2)
i i i i
Moreover,wepreparethefollowingadditionalconditionstoassureweakconsistencyandasymptotic
βˆ
normalityof ,whichcanbederivedbyslightlymodifyingtheresultsin[14]:
CC..54 ::lsiemquinefnmc→e{∞xi}λl(i(cid:2)esmi=in1Fniw(cid:3)iit(hβ()I⊗r ⊗xixx(cid:2)i)/(cid:2)βn)∈>(cid:5)00,,wfohrearlelλβ(A∈)Bin,dicatesthesmallesteigenvalueof
symmetricmatrixA.
βˆ
AccordingtoCorollary1in[14], hasw(cid:2)eakconsistencyandasymptoticnormalityunderthesecondi-
tions.Furthermore,fromConditionC.5, mi=1ni(cid:3)i(β)⊗xix(cid:2)i =O(n)issatisfied.Undertheassumption
βˆ
thatallconditionsaresatisfied, canbeformallyexpandedasfollows:
4332 H.Yanagiharaetal./LinearAlgebraanditsApplications436(2012)4329–4341
βˆ =β+ √1 b + 1b + √1 b +O (n−2), (3)
1 2 3 p
n n n n
whereb ,b ,andb arekr-dimensionalrandomvectors.Thepurposeofthissectionisachievedby
1 2 3
specifyingb ,b ,andb .
1 2 3
(cid:3)(β) β =βˆ (cid:3)(β)
Sincethelog-likelihoodfunction isamaximumat ,thefirstderivativeof becomes
β =βˆ
0 at ,i.e.,
kr
(cid:16)
∂(cid:3)(β)(cid:16)(cid:16)(cid:16) =(cid:5)m (cid:17)(y ⊗x)−n(p(βˆ)⊗x)(cid:18)=0 , (4)
∂β (cid:16) i i i i i kr
β=βˆ i=1
where0 isakr-dimensionalvectorofzeros.ToexpandEq.(4),wepreparethefollowingthreematrices
kr
−(cid:3)(β)/
consistingofthesecond,third,andfourthderivativesof n:
(cid:3) (cid:4)
1∂2(cid:3)(β) 1 ∂ ∂2
(β)=− , (β)=− ⊗ (cid:3)(β),
G2 ∂β∂β(cid:2) G3 ∂β(cid:2) ∂β∂β(cid:2)
n n
(cid:3) (cid:4)
1 ∂2 ∂2
(β)=− ⊗ (cid:3)(β).
G4 ∂β∂β(cid:2) ∂β∂β(cid:2)
n
(cid:3)(β) (β) (β)
Theresultofthefirstderivativeof in(4)impliesthefollowingexplicitformsofG ,G ,
2 3
(β)
andG (detailsofthederivationsaregiveninAppendixA):
4
(cid:19) (cid:20)
G2(β)=(cid:5)m ρi ∂p∂iβ(β(cid:2)) ⊗xi =(cid:5)m ρi(cid:21)(cid:3)i(β)⊗xix(cid:2)i(cid:22), (5)
i=1 i=1
(cid:19)(cid:3) (cid:4) (cid:20)
(cid:5)m ∂ ∂
(β)= ρ ⊗ (β) ⊗
G3 i ∂β(cid:2) ∂β(cid:2) pi xi
i=1
(6)
(cid:5)m (cid:17) (cid:18)
= ρi (cid:4)3,i(β)⊗xix(cid:2)i ,
i=1
(cid:19)(cid:3) (cid:4) (cid:20)
(cid:5)m ∂2 ∂
(β)= ρ ⊗ (β) ⊗
G4 i ∂β∂β(cid:2) ∂β(cid:2) pi xi
i=1
(7)
(cid:5)m (cid:17) (cid:18)
= ρi (cid:4)4,i(β)⊗xix(cid:2)i ,
i=1
where(cid:4)3,i(β)and(cid:4)4,i(β)arekr×(kr)2and(kr)2×(kr)2matrices,respectively,whicharedefined
by
(cid:5)r
(cid:4)3,i(β)= pia(β)e(cid:2)a⊗x(cid:2)i⊗qi,a(β)qi,a(β)(cid:2)−pi(β)(cid:2)⊗x(cid:2)i⊗(cid:3)i(β),
a=1
(cid:5)r
(cid:4)4,i(β)= pia(β)qi,a(β)qi,a(β)(cid:2)⊗xix(cid:2)i⊗(qi,a(β)qi,a(β)(cid:2)−pi(β)pi(β)(cid:2))
a=1 (8)
−(cid:3)(β)⊗xx(cid:2)⊗((cid:3)(β)−p(β)p(β)(cid:2))
i i i i i i
(cid:5)r
− (β) (β) (β) (β)(cid:2)⊗ (cid:2)⊗( (cid:2) + (cid:2)).
pia pib qi,a qi,b xixi eaeb ebea
a,b
H.Yanagiharaetal./LinearAlgebraanditsApplications436(2012)4329–4341 4333
(β)=
Here,ejisanr-dimensionaljth(cid:2)coordinateunit(cid:2)vectorwho(cid:2)sejthelementis1andothersare0,qi,a
ea−pi(β),andthenotation ra1,...,aj means ra1=1··· raj=1.
βˆ =β
ApplyingaTaylorexpansionaround toEq.(4)yields
1(cid:5)m {( − (β))⊗ } = (β)(βˆ −β)+ 1 (β)(cid:17)(βˆ −β)⊗(βˆ −β)(cid:18)
y np x G G
i i i i 2 3
n i=1 2
(cid:17) (cid:18) (cid:17) (cid:18)
+1 ⊗(βˆ −β)(cid:2) (β) (βˆ −β)⊗(βˆ −β)
I G
kr 4
6
+O (n−2). (9)
p
Noticethattheorderoftheleft-handsideofEq.(9)isO (n−1/2).BycomparingtheO (n−1/2),O (n−1),
p p p
andO (n−3/2)termsaftersubstituting(3)into(9),b ,b ,andb in(3)arespecifiedas
p 1 2 3
b = √1 G (β)−1(cid:5)m {(y −np(β))⊗x},
1 2 i i i i
n i=1
b =−1G (β)−1G (β)(b ⊗b ), (10)
2 2 3 1 1
2 (cid:23) (cid:24)
b =−1G (β)−1 G (β)(b ⊗b +b ⊗b )+1(I ⊗b(cid:2))G (β)(b ⊗b ) .
3 2 2 3 1 2 2 1 3 kr 1 4 1 1
βˆ
Weusethestochasticexpansionof withb ,b ,andb toevaluatethebiasoftheAICtotherisk
1 2 3
function.Thestochasticexpansionisregardedasaspecialcaseofthegeneralstochasticexpansionof
MLE,e.g.,in[15].
3. Mainresult
LetL(β)bealossfunctiondefinedby
L(β)=E[−2(cid:3)(β)]
⎧ ⎛ ⎞⎫
=− (cid:5)m ⎨( ∗⊗ )(cid:2)β− ⎝ +(cid:5)r ( (cid:2)β )⎠⎬, (11)
2 n p x log 1 exp x
i⎩ i i i j ⎭
i=1 j=1
∗
where p is the cell probability vector of the true model. Then, the risk function consisting of the
i
predictiveK–Linformationisgivenby
R=E[L(βˆ)]. (12)
Inthissection,weproposeaCAICthatimprovesthebiasoftheAICtoO(n−2)undertheassumption
thatthecandidatemodelincludesthetruemodel.NoticethatthecrudeAICisdefinedby
=− (cid:3)(βˆ)+ .
AIC 2 2kr (13)
− (cid:3)(βˆ)
Thus,itissufficienttoderivethebiasof 2 toRforevaluatingthebiasoftheAIC.Alsonotice
∗ = (β) − (cid:3)(βˆ)
thatp p holdswhenthecandidatemodelincludesthetruemodel.Then,thebiasof 2
i i
toRundertheassumptionthatthecandidatemodelincludesthetruemodelisexpandedas
= − [− (cid:3)(βˆ)]
B R E 2
(cid:5)m (cid:25) (cid:26)
= {( − (β))⊗ }(cid:2)βˆ
2 E y np x
i i i i
i=1
4334 H.Yanagiharaetal./LinearAlgebraanditsApplications436(2012)4329–4341
√ (cid:25) (cid:26)
= (cid:2) (β)βˆ
2 nE b G
(cid:23)√ 1 2
= [ (cid:2) (β)β]+ [ (cid:2) (β) ]
2 nE b G E b G b
1 2 1 2 1
(cid:20)
(14)
+√1 E[b(cid:2)G (β)b ]+ 1E[b(cid:2)G (β)b ] +O(n−2),
n 1 2 2 n 1 2 3
(β) (β) (β)
wherematricesG ,G ,andG aregivenby(5),(6),and(7),respectively,andkr-dimensional
2 3 4
randomvectorsb ,b ,andb aregivenby(10).Inmanycasesofpracticalinterest,amomentofstatistic
canbeexpanded1as2apowe3rseriesinn−1 (seee.g.[16,p.46]).Hence,theorderoftheremainder
termof(14)isshownbyO(n−2),notO(n−3/2).Indeed,ann−3/2termofthestochasticexpansionof
(cid:2)mi=1{(yi−nipi(β))⊗xi}(cid:2)βˆ inthebiascanbeexpressedasafifth-orderpolynomialofelements
ofb .Sinceb hasanasymptoticnormality,theexpectationofanodd-orderpolynomialofelements
ofb1 become1sO(n−1/2).Giventhisfact,theorderoftheremaindertermoftheexpansionin(14)is
O(n−12).
(cid:2) (β)
Fromelementarylinearalgebraandthedefinitionofb in(10),b G b in(14)isexpressedby
2 1 2 2
thefunctionofb as
1
(cid:2) (β) =−1 (cid:2) (β)( ⊗ )=−1 { (β)( ⊗ (cid:2))}.
b G b b G b b tr G b b b (15)
1 2 2 2 1 3 1 1 2 3 1 1 1
Sincethederivativeisinvarianttochangesintheorderofdifferentiation,wehave
(cid:2) (β)( ⊗ )= (cid:2) (β)( ⊗ )= (cid:2) (β)( ⊗ )=( ⊗ )(cid:2) (β)(cid:2) .
b G b b b G b b b G b b b b G b
1 3 1 2 1 3 2 1 2 3 1 1 1 1 3 2
Itfollowsfromtheaboveequationsandthedefinitionofb in(10)that
2
(cid:2) (β)( ⊗ + ⊗ )= ( ⊗ )(cid:2) (β)(cid:2)
b G b b b b 2 b b G b
1 3 1 2 2 1 1 1 3 2
=−(b ⊗b )(cid:2)G (β)(cid:2)G (β)−1G (β)(b ⊗b )
(cid:17)1 1 3 2 3 1 (cid:18)1
=−tr G (β)(cid:2)G (β)−1G (β)(b b(cid:2) ⊗b b(cid:2)) .
3 2 3 1 1 1 1
(cid:2) (β)
Thus,fromtheaboveresultandthedefinitionofb in(10),b G b in(14)isexpressedbythe
3 1 2 3
functionofb as
1
(cid:2) =−1 (cid:2) (β)( ⊗ + ⊗ )− 1 (cid:2)( ⊗ (cid:2)) (β)( ⊗ )
b G b b G b b b b b I b G b b
1 2 3 2 1 3 1 2 2 1 6 1 kr 1 4 1 1
(cid:17) (cid:18)
= 1tr G (β)(cid:2)G (β)−1G (β)(b b(cid:2) ⊗b b(cid:2)) − 1tr{G (β)(b b(cid:2) ⊗b b(cid:2))}. (16)
2 3 2 3 1 1 1 1 6 4 1 1 1 1
Hence,Eqs.(15)and(16)indicatethattheexpansionofBin(14)canbecalculateduntilthefourth
momentofb .
1
Sinceb consistsofacentralizedy,wecandirectlycalculatetheexpectationsin(14)bycentralized
1 i
,..., − (β)
momentsofy y .Then,allcombinationsofmultivariatemomentsofy np areneeded
1 m i i i
untilthefourth-order.However,itistroublesometocalculatethethird-andfourth-ordermultivariate
− (β)
momentsofy np ,becausewehavetoconsiderallcombinationsofthemultivariatemoments.
i i i
Forsimplicity,therelationsbetweenthemomentsofb andtheexpectationsofthederivativesof
1
−(cid:3)(β) β − (β)
withrespectto areusedinsteadofcalculatingthemultivariatemomentsofy np .It
i i i
[ ] = [ ] = (β) (cid:3)(β)
iseasytoobtainE b 0 becauseE y np .Fromtheresultofthefirstderivativeof
1 kr i i i
in(4)andthedefinitionofb in(10),wecanseethat
1
∂ √
(cid:3)(β)= (β) .
nG b
∂β 2 1
[∂(cid:3)(β)/∂β]= (β) (β) (β)
NoticethatE 0 holdsandG ,G ,andG areconstantmatrices.Byapplying
kr 2 3 4
generalformulasofexpectations(B.6)inAppendixBtothecaseofthemultinomiallogisticregression
H.Yanagiharaetal./LinearAlgebraanditsApplications436(2012)4329–4341 4335
model,thefollowingequationsareobtained:
(β)= (β) [ (cid:2)] (β),
nG nG E b b G
2 √2 1 1 2
(β)= (β) [ (cid:2) ⊗ (cid:2)]( (β)⊗ (β)),
nG n nG E b b b G G
3 2 1 1 1 2 2
nG (β)=n2(G (β)⊗G (β))E[b b(cid:2) ⊗b b(cid:2)](G (β)⊗G (β))
4 (cid:17)2 2 1 1 1 1 2 2 (cid:18)
−n2 (Ik2r2 +Kkr)(G2(β)⊗G2(β))+vec(G2(β))vec(G2(β))(cid:2) ,
( )
wherevec A isanoperatortotransformamatrixtoavectorbystackingthefirsttothelastcolumn
( )=( (cid:2),..., (cid:2) )(cid:2) =( ,..., )
ofA,i.e.,vec A a a whenA a a (seee.g.[17,Chapter16.2]),andK isthe
1 m 1 m m
m2×m2vec-permutationmatrixsuchthatvec(B)=K vec(B(cid:2))whenBisanm×mmatrix(seee.g.
m
[17,Chapter16.3]).Theseresultsleadustothesimpleexpressionofmomentsofb as
1
E[b b(cid:2)]=G (β)−1, (17)
1 1 2
E[b(cid:2) ⊗b b(cid:2)]= √1 G (β)−1G (β)(G (β)−1⊗G (β)−1), (18)
1 1 1 n 2 3 2 2
E[b1b(cid:2)1⊗b1b(cid:2)1]=+(Ikv2erc2(+G K(βkr))−(1G)2v(eβc)(−G1(⊗β)G−21()β(cid:2))+−1O)(n−1). (19)
2 2
Theresultin(17)impliesthat
E[b(cid:2)G (β)b ]=E[{G (β)b b(cid:2)}]=tr{G (β)G (β)−1}=kr. (20)
1 2 1 2 1 1 2 2
Similarly,from(15)and(18),wehave
[ (cid:2) (β) ]=−1 [ { (β)(cid:2)( (cid:2) ⊗ (cid:2))}]
E b G b E tr G b b b
1 2 2 2 3 1 1 1
(cid:17) (cid:18)
=− √1 tr G (β)(cid:2)G (β)−1G (β)(G (β)−1⊗G (β)−1) . (21)
3 2 3 2 2
2 n
(β) = (β)
NoticethatG K G holdsbecausethederivativeisinvarianttochangesintheorderof
3 kr 3
differentiation.ByusingthisfactandEq.(19),theexpectationofthefirstpartin(16)isgivenby
(cid:25) (cid:17) (cid:18)(cid:26)
E tr G (β)(cid:2)G (β)−1G (β)(b b(cid:2) ⊗b b(cid:2))
3 2 3 1 1 1 1
(cid:17) (cid:18)
=tr G3(β)(cid:2)G2(β)−1G3(β)(Ik2r2 +Kkr)(G2(β)−1⊗G2(β)−1)
+vec(G (β)−1)(cid:2)G (β)(cid:2)G (β)−1G (β)vec(G (β)−1)+O(n−1)
(cid:17) 2 3 2 3 2 (cid:18)
=2tr G (β)(cid:2)G (β)−1G (β)(G (β)−1⊗G (β)−1)
3 2 3 2 2
(22)
+vec(G (β)−1)(cid:2)G (β)(cid:2)G (β)−1G (β)vec(G (β)−1)+O(n−1).
2 3 2 3 2
Moreover,sincethederivativeisinvarianttochangesintheorderofdifferentiation,wecanseethat
(β) = (β)
G K G and
4 kr 4
(cid:17) (cid:18)
tr G (β)(G (β)−1⊗G (β)−1) =vec(G (β)−1)(cid:2)G (β)vec(G (β)−1).
4 2 2 2 4 2
ByusingtheaboverelationsandEq.(19),theexpectationofthesecondpartin(16)isgivenby
(cid:25) (cid:26)
E tr{G4(β)(b1b(cid:2)1⊗b1b(cid:2)1)} =tr{G4(β)(Ik2r2 +Kkr)(G2(β)−1⊗G2(β)−1)}
+vec(G (β)−1)G (β)vec(G (β)−1)+O(n−1)
2 4 2
=3tr{G (β)(G (β)−1⊗G (β)−1)}+O(n−1). (23)
4 2 2
4336 H.Yanagiharaetal./LinearAlgebraanditsApplications436(2012)4329–4341
Hence,fromEqs.(16),(22),and(23),wecanseethat
(cid:25) (cid:26) (cid:17) (cid:18)
E b(cid:2)G b =tr G (β)(cid:2)G (β)−1G (β)(G (β)−1⊗G (β)−1)
1 2 3 3 2 3 2 2
+ 1vec(G2(β)−1)(cid:2)G3(β)(cid:2)G2(β)−1G3(β)vec(G2(β)−1) (24)
2
(cid:17) (cid:18)
− 1tr G (β)(G (β)−1⊗G (β)−1) +O(n−1).
4 2 2
2
[ (cid:2) (β)β]=
Consequently,bysubstitutingE b G 0,andEqs.(20),(21),and(24)into(14),thebiasof
1 2
− (cid:3)(βˆ)
2 toRisexpandedas
B=2kr+ 1{α (β)+α (β)−α (β)}+O(n−2),
1 2 3
n
α (β) α (β) α (β)
wherecoefficients , ,and aregivenby
1 2 3
(cid:17) (cid:18)
α (β)=tr G (β)(cid:2)G (β)−1G (β)(G (β)−1⊗G (β)−1) ,
1 3 2 3 2 2
α (β)=vec(G (β)−1)(cid:2)G (β)(cid:2)G (β)−1G (β)vec(G (β)−1), (25)
2 (cid:17) 2 3 2 3(cid:18) 2
α (β)=tr G (β)(G (β)−1⊗G (β)−1) .
3 4 2 2
− (cid:3)(βˆ)
TheCAICcanthenbedefinedbyaddinganestimatedBto 2 ,i.e.,
(cid:17) (cid:18)
=− (cid:3)(βˆ)+ + 1 α (βˆ)+α (βˆ)−α (βˆ) .
CAIC 2 2kr (26)
1 2 3
n
TheCAICimprovesthebiasoftheAICtoO(n−2),althoughtheorderofthebiasoftheAICisO(n−1),
i.e.,thefollowingequationsaresatisfied:
R−E[AIC]=O(n−1), R−E[CAIC]=O(n−2),
whereRistheriskfunctiongivenby(12).
4. Numericalstudies
Inthissection,weconductnumericalstudiestoshowthattheCAICin(26)worksbetterthanthe
crudeAICin(13).TocomparetheperformancesoftheAICandtheCAIC,thefollowingtwoproperties
areconsidered:
(I) theselectionprobability:thefrequencyofthemodelchosenbyminimizingtheinformation
criterion.
(II) thepredictionerrorofthebestmodel(PE ):theriskfunctionofthebestmodelchosenbythe
B
informationcriterion,whichisdefinedby
PE =E[L(βˆ )],
B B
whereL(β)isthelossfunctiongivenby(11)andβˆ istheMLEofβunderthebestmodel.
B
ThesetwopropertieswereevaluatedbyaMonteCarlosimulationwith10,000iterations.Theinforma-
tioncriterionwiththehigherselectionprobabilityofthetruemodelandthesmallerpredictionerror
ofthebestmodelisregardedasahigh-performancemodelselector.InthebasicconceptoftheAIC,a
goodmodelselectionmethodisonethatchoosesthebestmodelsothatthepredictionisimproved.
Hence,PE isamoreimportantpropertythantheselectionprobability.
B
,..., = = ( = ,..., )
WepreparedeightcandidatemodelsM M ,withm 20and50,n 5 i 1 m and
= × 1 =8( ,..., )(cid:2) i
r 2.Anm 8matrixofexplanatoryvariablesX x x wasconstructedasfollows.Thefirst
1 m
H.Yanagiharaetal./LinearAlgebraanditsApplications436(2012)4329–4341 4337
Table1
Selectionprobabilityofthemodelandthepredictionerrorofthebestmodel.
Casem Criterion Selectionprobability PEB
M1 M2 M3 M4 M5 M6 M7 M8
1 20 AIC 1.81 0.29 74.84 11.41 4.95 2.88 2.20 1.62 210.06
CAIC 3.19 0.66 79.92 10.01 3.58 1.38 0.88 0.38 207.84
50 AIC 0.01 0.00 79.15 11.16 4.77 2.21 1.56 1.14 511.32
CAIC 0.01 0.00 81.25 10.71 4.27 1.88 1.17 0.71 511.04
2 20 AIC 77.22 10.92 4.86 2.69 1.56 1.11 0.77 0.87 202.42
CAIC 82.63 10.06 3.81 2.07 0.76 0.38 0.19 0.10 200.41
50 AIC 79.21 10.89 4.48 2.20 1.24 1.04 0.55 0.39 494.89
CAIC 80.99 10.58 4.10 1.91 1.01 0.70 0.44 0.27 494.63
Note:Theselectionprobabilityofthetruemodelismarkedinbold.
columnofXis1 ,where1 isanm-dimensionalvectorofones,andtheremainingsevencolumnsof
m m
( , . )
XweregeneratedrandomlyfromthebinomialdistributionB 1 05 .Simulationdataweregenerated
β∗ = (β∗(cid:2),β ∗(cid:2))(cid:2)
fromthemultinomialdistributionwiththetruecellprobabilityconsistingof .In
β∗ 1 2
thissimulationstudy,wepreparedtwo ,asfollows:
: β∗ =( , . ,− . , , , , , )(cid:2), β∗ =(− . ,− . , . , , , , , )(cid:2),
Case1 0 02 10 0 0 0 0 0 01 04 12 0 0 0 0 0
1 2
: β∗ =(− . , , , , , , , )(cid:2), β∗ =( . , , , , , , , )(cid:2).
Case2 05 0 0 0 0 0 0 0 07 0 0 0 0 0 0 0
1 2
( = ,..., )
ThematrixofexplanatoryvariablesinM consistsofthefirstjcolumnsofX j 1 8 .Thus,the
j
truemodelinCase1isM ,andthetruemodelinCase2isM .
3 1
Table1showsthetwoproperties(I)and(II).Inthetable,theselectionprobabilityofthetruemodel
ismarkedinbold.Fromthistable,wecanseethattheselectionprobabilitiesandthepredictionerrors
oftheCAICwereimprovedincomparisonwiththoseoftheAICinallsituations.Wesimulatedseveral
othermodelsandobtainedsimilarresults.
5. Conclusionanddiscussion
Inthispaper,weproposedtheCAICforselectingvariablesinthemultinomiallogisticregression
models.TheproposedCAICimprovesthebiasoftheAICtoO(n−2),althoughtheorderofthebiasof
theAICisO(n−1).Byusingrelationsbetweenthemomentsofb andexpectationsofthederivatives
1
−(cid:3)(β)
of insteadofdirectlycalculatingthemomentsofy toevaluatethemomentsofb ,asimple
i 1
expressionoftheCAICisdeveloped.Indeed,thebiascorrectiontermoftheproposedCAICisrepre-
(βˆ) (βˆ) (βˆ)
sentedbyonlythreematricesG ,G ,andG ,whichconsistofthesecond,third,andfourth
2 3 4
−(cid:3)(β) (βˆ) (βˆ) (βˆ)
derivativesof .EventhoughexpressionsofG ,G ,andG arenotsimple,wecan
2 3 4
derivethebiascorrectiontermoftheCAICfromlinearfunctionsofG (βˆ)−1,G (βˆ),andG (βˆ).This
2 3 4
isadesirablecharacteroftheCAIC.
Inallsituationsofthesimulationstudy,theCAICimprovedthecrudeAICinthesenseofmaking
ahighselectionprobabilityofthetruemodelandasmallpredictionerrorofthebestmodelchosen
bytheinformationcriterion.However,theimprovementsweresmallerwhenthesamplesizewas
large.ThisisnaturalbecausetheCAICisproposedsothatthebiasoftheAICiscorrectedwhenthe
samplesizeissmall.Needlesstosay,theAICandtheCAICareasymptoticalequivalents.Hence,the
differencebetweentwocriteriabecomessmallwhenthesamplesizeisincreased.Thesamplesizes
ofoursimulationwere100and250.Nevertheless,acleardifferenceexistsintheperformancesofthe
CAICandtheAIC.ThisdifferenceindicatesthattheCAICisvaluableevenwhenthesamplesizeisnot
sosmall.Consequently,werecommendusingtheCAICinsteadoftheAICforselectingmultinomial
logisticregressionmodels.
ThesimpleexpressionoftheproposedCAICisbasedonthepropertythatthesecondderivatives
−(cid:3)(β)
of donotdependonresponsevariables.Ageneralizedlinearmodel(GLM)withanaturallink
andaknowndispersionparameter,e.g.,alogisticregressionmodeloraPoissonregressionmodel,will
4338 H.Yanagiharaetal./LinearAlgebraanditsApplications436(2012)4329–4341
havethisproperty.Then,wecansimplyexpressthebias-correctedAICjustliketheproposedCAIC
in(26)inthesamewaypresentedinSection3.Namely,thebias-correctedAICwithconstantsecond
derivativesofthenegativelog-likelihoodfunctionmaybestatedby
= +γ (θˆ)+γ (θˆ)−γ (θˆ),
CAIC AIC
1 2 3
θˆ θ γ (θ) γ (θ) γ (θ)
where istheMLEofunknownparameter ,andcoefficients , ,and aregivenby
1 2 3
(cid:17) (cid:18)
γ (θ)=tr C(θ)(cid:2)H(θ)−1C(θ)(H(θ)−1⊗H(θ)−1) ,
1
γ (θ)=vec(H(θ)−1)(cid:2)C(θ)(cid:2)H(θ)−1C(θ)vec(H(θ)−1),
2 (cid:17) (cid:18)
γ (θ)=tr Q(θ)(H(θ)−1⊗H(θ)−1) .
3
(θ) (θ) (θ)
Here,H ,C ,andQ arematricesconsistingofthesecond,third,andfourthderivatives,respec-
tively,ofthenegativelog-likelihoodfunctionandaredefinedby(B.4)inAppendixB.
Beforeconcludingthissection,weconsideralog-likelihoodratiostatisticfortestingforthenull
:β =β
hypothesisH .Itisknownthatalog-likelihoodratiostatisticfortestingthenullhypothesis
0 0
= {(cid:3)(βˆ)−(cid:3)(β )}
H isT 2 .BycomparingEq.(25)withthegeneralformulaoftheBartlettfactorofa
0 0
log-likelihoodratiostatistic(seee.g.[15]),wefindthatthefirsttermintheasymptoticexpansionof
thebiasofAICconsistsofthesimilarcoefficientsinaBartlettfactorofT.Thisisbecausethebiasof
[(cid:3)(βˆ)] = [(cid:3)(βˆ)]+
AICandtheBartlettfactorofT arepartiallyformedfrom2E .However,B 2E Rand
[ ]= [(cid:3)(βˆ)]− [(cid:3)(β )] − [(cid:3)(β )]
E T 2E 2E arecleardifferentbecauseRisnotequalto 2E .Althoughboth
0 0
[ ]
aredifferent,wewillobtainanasymptoticexpansionofE T inthesamewayinSection3.Then,the
BartlettfactorofT maybeexpressedbythelinearfunctionsofG (βˆ)−1,G (βˆ)andG (βˆ)asinthe
2 3 4
biascorrectiontermofCAIC.
(β) (β) (β)
AppendixA. ExplicitformsofG ,G ,ANDG
2 3 4
Inthissubsection,forsimplicity,wewrite(cid:3)(β),p(β),andp (β)as(cid:3),p,andp ,respectively.
i i ij i i ij
Noticethat
∂pi = (cid:27) − (cid:28) (cid:2), ( = ,..., ),
∂β(cid:2) pij ej pi xi j 1 r
j
wheree isthejthcoordinateunitvector,whichisusedinEq.(8).ThisresultandEq.(2)implythat
j
∂ (cid:21) (cid:22)
∂βpi(cid:2) = pi1(e1−pi)x(cid:2)i,...,pir(er −pi)x(cid:2)i =(cid:3)i⊗x(cid:2)i.
(β)
SubstitutingtheaboveresultintothedefinitionofG yieldsEq.(5).Furthermore,fromthedefini-
2
tionsofG3(β)andG4(β),wecanseethat(cid:4)3,i(β)and(cid:4)4,i(β)in(6)and(7),respectively,satisfy
∂ ∂2
(cid:4) (β)= ⊗(cid:3), (cid:4) (β)= ⊗(cid:3).
3,i ∂β(cid:2) i 4,i ∂β∂β(cid:2) i
Noticethatthe(a,b)thelementof(cid:3) isp δ −p p ,whereδ istheKroneckerdelta,i.e.,δ =1
i ia ab ia ib ab aa
andδab =0fora(cid:11)=b.Thisequationleadsustootherexpressionsof(cid:4)3,i(β)and(cid:4)4,i(β),asfollows:
(cid:5)r ∂
(cid:4)3,i(β)= ∂β(cid:2)(piaδab−piapib)⊗eae(cid:2)b,
a,b
(A.1)
(cid:5)r ∂2
(cid:4)4,i(β)= ∂β∂β(cid:2)(piaδab−piapib)⊗eae(cid:2)b.
a,b