Some remarks on Bayesian inference for one-way ANOVA models

16 Pages·2008·0.19 MB·English

AnnInstStatMath(2008)60:483–498
DOI10.1007/s10463-007-0117-5
Some remarks on Bayesian inference for one-way
ANOVA models
FabrizioSolari · BruneroLiseo · DongchuSun
Received:14July2005/Revised:13November2006/Publishedonline:11May2007
©TheInstituteofStatisticalMathematics,Tokyo2007
Abstract Weconsiderthestandardone-wayANOVAmodel;itiswell-known
thatclassicalstatisticalproceduresarebasedonascalarnon-centralityparame-
ter.Inthispaperweexplorebothmarginallikelihoodandintegratedlikelihood
functions for this parameter and we show that they exactly lead to the same
answer. On the other hand, we prove that a fully Bayesian testing procedure
mayprovidedifferentconclusions,dependingonwhatisconsideredtobethe
realquantityofinterestinthemodelor,saiddifferently,whicharethecompet-
inghypotheses.Weillustratetheseissuesviaarealdataexample.
Keywords Integrated likelihood · Marginal likelihood · Model choice ·
ObjectiveBayesfactor·Referenceprior
1 Introduction
Analysis of variance (ANOVA) is an extremely important method in explor-
atory and confirmatory data analysis (Gelman, 2005). Here we focus on the
one-way ANOVA, where we assume that data are observed according to the
(full)modelM
F
F.Solari
ISTAT,ViaMagenta4,Rome00184,Italy
B
B.Liseo( )
Dip.distudigeoeconomici,linguistici,statisticiestorici,
UniversitàdiRomaSapienza,vialedelCastroLaurenziano9,Rome00161,Italy
e-mail:[email protected]
D.Sun
DepartmentofStatistics,UniversityofMissouri,
146MiddlebushHall,Columbia,MO65211-6100,USA
484 F.Solarietal.
Y =µ +(cid:1) , j=1,...,n, i=1,...,k, (1.1)
ij i ij i
where µ = (µ ,...,µ ) and σ2 are unknown parameters and (cid:1) are i.i.d.
1 k ij
N(0,σ2).Classicalanalysisofvariancefocusesonthehypothesestest
M :µ =···=µ versus M :µ (cid:1)=µ foratleastonepair (i,j) (1.2)
0 1 k 1 i j
anditisusuallybasedonthestatistic
(cid:1)
F = (cid:1) ki=(cid:1)1ni(Y¯i−Y¯)2/(k−1) ,
ki=1 nj=i1(Yij−Y¯i)2/(n−k)
(cid:1) (cid:1) (cid:1) (cid:1)
wheren= ki=1ni,Y¯i =n−i 1 nj=i1Yij andY¯ =n−1 ki=1 nj=i1Yij.
It is well-known that under the null hypothesis M , the F statistic follows
0
a Fk−1,n−k distribution. Under the alternative hypothesis M1, the F statistic
is distributed according to a non-central Fk−1,n−k,λ distribution, whose non-
centralityparameterλisgivenby
(cid:2)k
1
λ= n(µ −µ¯)2,
σ2 i i
i=1
(cid:1)
withµ¯ =n−1 ki=1niµi.
The above testing procedure is a clear and classical example of how the
problemofeliminatingnuisanceparametersishandledinfrequentiststatistics.
AlthoughthefullmodelM hask+1parameters,ascalarteststatisticiscon-
F
structedtocomparedifferentvaluesofthescalarnon-centralityparameterλ.In
acertainsense,theclassicaltestactsasifwewouldhaveobservedonlythemar-
ginalexperiment,whichprovidesthequantityF,discardingalltheinformation
aboutthesinglemeantreatments,(y¯ ,...,y¯ ).
1 k
Itisthenapparentthatthereisalossofinformationinthisprocedure.How-
ever,itisnotclearhowmuchinformationislost,andthispaperisanattemptto
clarifytheseissues.Bertolinoetal.(1990)haveshownthattheclassicaltestcan
bereinterpretedintermsofamarginallikelihoodoftheparameterτ =λ/n;we
willbrieflyrecalltheirapproachinSect.2.InSect.3weperformanintegrated
likelihoodapproachtotheproblemandshowthat,usingaconditionalreference
prior(BergerandBernardo1992),theresultingintegratedlikelihoodisexactly
equal to the marginal likelihood of Bertolino et al. (1990). This equivalence
amountstosay,inouropinion,thatwhenτ istherealquantityofinterest,the
classicaltest,oratleastitsmarginallikelihoodcounterpart,isacorrectreport
oftheinformationavailable.InSect.4wederivetheclassoffirst-ordermatch-
ingpriorsfortheparameterofinterestandperformsomecomparisonsamong
different priors in terms of frequentist coverage. In Sect. 5 we also perform
anobjectiveBayesanalysisforcomparingdifferentmodels.Weshowthatthe
classicalF test(anditslikelihoodcounterpart)correspondtoaspecificmodel
Bayesianone-wayANOVA 485
comparisonthatmaybeormaybenottheappropriatetestinone-wayANOVA
problems;weillustrateourfindingsviaasimplerealdataexample.
2 Marginallikelihoodapproach
ThegreatsimplicityofthestandardANOVAmodelreliesonthetransforma-
tion of a multiparametric problem into an equivalent scalar one, where the
only parameter involved is the non-centrality parameter λ, i.e. rewriting the
hypothesestest(1.2)as
M :λ=0 versus M :λ>0. (2.1)
0 1
Then, frequency-based solutions to the ANOVA problem consist only in the
computationofthesamplerealizationF oftheFstatistic.Ontheotherhand,
obs
Bayesiansolutionsarecomplicatedbecauseofthepossiblyhighdimensionof
the parameter (µ,σ2). Bertolino et al. (1990) propose to consider F as the
obs
actualresultofamarginalexperimentandthentousethesamplingdistribution
oftheF statisticasamarginallikelihood.
Theyalsosuggestthatitismoreconvenienttoworkwiththeparameter
λ
τ =
n
insteadofλ.Thischoicehastwoadvantages.First,τ hasaclearinterpretation
astheratioofbetweenandwithinvariances;second,itdependsonthesample
sizenonlythroughtherelativefrequenciesn/n.Denotingbyd = k−1and
i 1
d =n−kthedegreesoffreedomoftheFstatistic,Bertolinoetal.(1990)write
2
themarginallikelihoodofτ as
(cid:3) (cid:4)
L(τ)∝(cid:2)∞ p (cid:3)nτ(cid:4)(cid:5) d(cid:3)1+2d2 +(cid:4)j (cid:5) d1Fobs (cid:6)j, (2.2)
j=0 j 2 (cid:5) d21 +j d2+d1Fobs
where,forj=1,2,...,p(z)=exp{−z}zj/j!.
j
Itisuseful,forreasonsthatwillbeclearlater,torewriteL(τ)asafunctionof
the Kummer confluent Hypergeometric function M(a,b,z) (Abramowitz and
Stegun1964);from(2.2)itimmediatelyfollowsthat
(cid:2)∞ (cid:3)nτ(cid:4)(cid:5)d +d (cid:6) (cid:5)d (cid:6)−1(cid:5) d F (cid:6)j
L(τ)∝ p 1 2 2 1 obs
j 2 2 2 d +d F
j=0 j j 2 1 obs
(cid:7) (cid:8) (cid:5) (cid:6)
nτ d +d d nτ d F
∝exp − M 1 2, 2, 1 obs ,
2 2 2 2 d +d F
2 1 obs
where(a) =(cid:5)(a+j)/(cid:5)(a)isthePochhammer’ssymbol.
j
486 F.Solarietal.
Making use of an asymptotic approximation for L(τ), Solari (2002) sho√ws
that a noninformative prior for the parameter τ is given by π(τ) ∝ 1/ τ.
Bertolino and Racugno (1994) perform a robust Bayesian analysis based on
the marginal likelihood L(τ). Their results seem to be strongly dependent on
thechoiceofthefunctionoftheparameter(µ,σ2)usedforreformulatingthe
hypothesestest(1.2)intheform(2.1)(Solari2002).
3 Integratedlikelihoodapproach
Inthissectionwederiveanintegratedlikelihoodfunctionfortheparameterτ.
Todoso,inSect.3.1wefirstproposeaone-to-onereparameterizationof(µ,σ2),
namely (τ,ξ), where ξ is a vector of suitable nuisance parameters; second, in
Sect.3.2,weapplythereferenceprioralgorithm.Weshowthat,nomatterwhat
istheorderofimportanceoftheparameters,themarginalreferencepriorfor
thenuisanceparameterξ isalwaysthesame,whilethemarginalprioronτ may
change;finallyinSect.3.3weobtaintheintegratedlikelihoodfunction
(cid:9)
L˜(τ)∝ L(τ,ξ)π(ξ|τ)dξ
ξ
inaclosedform.
3.1 Reparameterizationof(µ,σ2)
Firstwere-expressmodel(1.1)inmatrixnotation.Let1 and0 denote,respec-
p p
tively,thep×1vectorsof1’sand0’s.Thus,
y=Aµ+(cid:2),
wherethen×kdesignmatrixAisgivenby
⎡ ⎤
1 0 ··· 0
⎢ n1 n1 n1⎥
A=⎣ .. .. .. .. ⎦
. . . .
0 0 ··· 1
nk nk nk
and(cid:2) isN(0,σ2I ).Then,theparameterτ canbere-writtenas
n
1 1
τ = (µ−µ¯1 )tR(µ−µ¯1 )= µtDtRDµ,
σ2 k k σ2
whereR=n−1 diag(n ,n ,...,n )and
1 2 k
Bayesianone-wayANOVA 487
⎡ ⎤
n−n −n ··· −n
1 2 k
1 ⎢⎢ −n1 n−n2 ··· −nk ⎥⎥
D= ⎢ . . . . ⎥.
n⎣ .. .. .. .. ⎦
−n −n ··· n−n
1 2 k
Note that the parameter of interest can be also expressed as τ = (nσ2)−1
(cid:1)
ki=1niµi(µi−µ¯)or,usingmatrixnotation,
1 1
τ = µtR(µ−µ¯1 )= µtRDµ.
σ2 k σ2
Nowconsiderthefollowingsteps.
Step1Sincethek-dimensionalvector(µ−µ¯1 )liesina(k−1)-dimensional
k
subspace of Rk, we need to transform (µ−µ¯1 ) into a (k−1)-dimensional
k
vector θ = (θ1,...,θk−1)t. As rank(RD) = rank(D) = k−1, there exists a
(k−1)×kmatrixPsuchthatPtP =RD.Settingθ =Pµ,wehave
Ptθ =PtPµ=RDµ=R(µ−µ¯1 )
k
andhence
µ=R−1Ptθ +µ¯1 .
k
TheJacobianmatrixH ofthetransformationfrom(µ,σ2)to(θ,σ2,µ¯)isthen
1
⎡ ⎤
H =(cid:16) ∂(µ,σ2) (cid:17)=⎣PR0t−1 0k1−1⎦.
1 ∂(θ,σ2,µ¯) k
1t 0
k
Step2Now,switchingtopolarcoordinatesasinBergeretal.(1998),wecan
(cid:1)
reparameterize θ in the new parameters η = θtθ = n−1 ki=1ni(µi −µ¯)2 and
ψ =(ψ1,...,ψk−2)∈(cid:5) ≡(0,π)k−3×(0,2π),thatis,
√
θ = ηcosψ ,
1 √ 1
θ = ηsinψ cosψ ,
2 1 2
.
.
.
√
θk−2 =√ηsinψ1···sinψk−3cosψk−2,
θk−1 = ηsinψ1···sinψk−3sinψk−2.
488 F.Solarietal.
Then,wehave
⎡ ⎤
(cid:16) (cid:17) η−1/2bt 0 0
H2 = ∂(∂η(,θσ,σ2,2µ,¯µ,¯ψ)) =⎢⎢⎣ 00tktk−−11 10 01 ⎥⎥⎦,
η1/2C 0k−2 0k−2
wherethe(k−1)-dimensionalvectorbandthe(k−2)×(k−1)matrixCare
givenby
b=η1/2 ∂∂θη = 12(cosψ1,...,sinψ1...sinψk−2)t
and
⎡ ⎤
∂θ ⎢−sinψ1 ··· cosψ1sinψ2×···×sinψk−2 ⎥
C=η−1/2 =⎣ .. .. .. ⎦.
∂ψ . . .
0 ··· sinψ1×···×sinψk−3cosψk−2
Step 3 Since τ = η/σ2, the Jacobian matrix H of the transformation from
3
(η,σ2,µ¯,ψ)to(τ,σ2,µ¯,ψ)isgivenby
⎡ ⎤
H3 =(cid:16)∂∂((τη,,σσ22,,µµ¯¯,,ψψ))(cid:17)=⎢⎢⎣ στ02 100 001 000tktktk−−−222⎥⎥⎦.
0k−2 0k−2 0k−2 Ik−2
We have defined a new set of parameters (τ,σ2,µ¯,ψ), in which τ is the
parameter of interest and ξ = (σ2,µ¯,ψ) is the nuisance parameter. It must
be noticed that the ordering in ξ is just one among the possible ordering we
couldhavechosen.However,sinceµ¯ andψcontaininformationaboutlocation
and direction of µ respectively, they do not seem to give direct information
abouttheparameterofinterestτ.Ontheotherhand,thereisastrongrelation
betweenτ andσ2.Consequently,itseemsreasonabletoconsidertheparameter
ordering{τ,σ2,µ¯,ψ},asweactuallydid,or,alternatively,{τ,σ2,ψ,µ¯}.Inany
case,itwillbeshownthatthereferencepriorsdonotdependontheordering
ofthecomponentsofthenuisanceparameterξ.
Forthenewparameter(τ,σ2,µ¯,ψ),theinformationmatrixI(τ,σ2,µ¯,ψ)is
givenby
I(τ,σ2,µ¯,ψ)=H H H I(µ,σ2)HtHtHt,
3 2 1 1 2 3
Bayesianone-wayANOVA 489
where
(cid:18) (cid:19) (cid:18) (cid:19)
1 AtA 0 1 R 0
I(µ,σ2)= σ2 k =n σ2 k .
0t n 0t 1
k 2(σ2)2 k 2(σ2)2
Somealgebrayields
(cid:20)(cid:18) (cid:19) (cid:21)
btb btb
I(τ,σ2,µ¯,ψ)=ndiag τ σ2 , 1 , τCCt ,
btb τbtb+1/2 σ2
σ2 (σ2)2
where CCt = diag(1,sin2ψ1,...,sin2ψ1···sin2ψk−3). Noting that btb = 1/4,
weget
(cid:20)(cid:18) (cid:19) (cid:21)
I(τ,σ2,µ¯,ψ)=ndiag 411τ τ4+σ122 , σ12, τCCt (3.1)
(cid:5)(cid:16)4σ2 4(σ(cid:17)2)2 (cid:6)
i i
≡ndiag 11 12 , i , I .
i i 33 44
21 22
where,thefirstelementofthediagonalinthelastformularefersto(τ,σ2),the
secondelementreferstoµ¯ andI referstoψ.
44
3.2 Referencepriors
3.2.1 Directapproach
Herewederivetheone-at-a-timereferencepriorfor(τ,σ2,µ¯,ψ)intheorder-
ing{τ,σ2,µ¯,ψ},inthesensethatτ istheparameterofinterestand(σ2,µ¯,ψ)
arethenuisanceparametersindescendingorderofinterest.Followingtheref-
erenceprioralgorithm(BergerandBernardo,1992),theconditionalreference
priorsfor(τ,σ2,µ¯,ψ)aregivenby
(cid:22) (cid:22) k(cid:23)−3
π (ψ|τ,σ2,µ¯)∝(cid:22)I (cid:22)12 ∝ (sinψ)k−i−2,
d 44 i
i=1
(cid:22) (cid:22)
π (µ¯|τ,σ2)∝(cid:22)i (cid:22)12 ∝ 1,
d 33
(cid:22) (cid:22)
π (σ2|τ)∝(cid:22)i (cid:22)12 ∝ 1 ,
d 22 σ2
(cid:22) (cid:22) (cid:22) (cid:22)
π (τ)∝(cid:22)i (cid:22)−12(cid:22)i i −i i (cid:22)12 ∝ √ 1
d 22 11 22 12 21 τ(τ +2)
anditiseasytocheckthattheydonotchangewhateveristheorderingofthe
nuisanceparameters(σ2,µ¯,ψ).
490 F.Solarietal.
Inalltheabovecomputations,wehaveconsideredthevectorparameterψas
asingleblock.Ifwesplititintoitsk−2distinctcomponents,thenewreference
priorforeachψ,i=1,...,k−2,hasthefollowinginterestingform:
i
(cid:23)i−1
πd(ψi |τ,σ2,µ¯,ψ1,...,ψi−1)∝ sinψj ∝ 1.
j=1
3.2.2 Reverseapproach
Now,wewillconstructthereversereferenceprior,wherethereverseordering
of the parameter, namely {ψ,µ¯,σ2,τ} is considered in the derivation of the
priordistributions.Thisway,thepriordistributionoftheparameterofinterest
τ conditionallyonthenuisanceparameter(ψ,µ¯,σ2)willbecomputedfirst.In
thiscase,wehave
(cid:5) (cid:16) (cid:17)(cid:6)
i i
I(ψ,µ¯,σ2,τ)=ndiag I , i , 11 12
44 33 i i
(cid:20) 2(cid:18)1 22 (cid:19)(cid:21)
1 τ+2 1
=ndiag τCCt, , 4(σ2)2 4σ2 .
σ2 1 1
4σ2 4τ
Itfollowsthattheconditionalreferencepriorsare
(cid:22) (cid:22)
π (τ |ψ,µ¯,σ2)∝(cid:22)i (cid:22)12 ∝ √1 ,
r 11 τ
(cid:22) (cid:22)
π (σ2 |ψ,µ¯)∝(cid:22)i (cid:22)12 ∝ 1 ,
r 22 σ2
(cid:22) (cid:22)
π (µ¯ |ψ)∝(cid:22)i (cid:22)12 ∝ 1,
r 33
(cid:22) (cid:22) k(cid:23)−3
π (ψ)∝(cid:22)I (cid:22)12 ∝ (sinψ)k−i−2.
r 44 i
i=1
Asbefore,thefinalanswerisindependentoftheorderingofthenuisancepara√m-
eters.Notethatthereverseapproachleadstothemarginalpriorπ (τ)∝1/ τ,
r
i.e.thesamepriorderivedinSolari(2002)forthemarginallikelihood(2.2).
Again,ifoneconsidersk−2groupsofparametersψ,i=1,...,k−2,instead
i
ofthesinglegroupψ,aconstantreferencepriorfortheψ’sisobtained
i
(cid:23)i−1
πr(ψi |ψ1,...,ψi−1)∝ sinψj ∝ 1.
j=1
Bayesianone-wayANOVA 491
3.3 Integratedlikelihood
In this section we will compute the integrated likelihood using the reference
priorsobtainedinSects.3.2.1and3.2.2.Itmustbenoticedthatboththedirect
andthereversereferencepriorproducethesamepriorforthenuisanceparam-
eter,i.e.
k(cid:23)−3
π(σ2,µ¯,ψ |τ)=π(σ2,µ¯,ψ)∝ 1 (sinψ)k−i−2,
σ2 i
i=1
resultinginthesameintegratedlikelihood
(cid:9)
L˜(τ)= L(τ,σ2,µ¯,ψ)π(σ2,µ¯,ψ)dσ2dµ¯ dψ.
Startingfromthelikelihoodfunctionfor(µ,σ2)
(cid:24) (cid:25)
1 1
L(µ,σ2)∝ exp − (y−Aµ)t(y−Aµ)
(σ2)n2 2σ2
andpassingthroughthereparameterizationfromStep1toStep3,aftersome
algebraweobtain
(cid:24)
1 nτ n(µ¯ −y¯)2 SST
L(τ,σ2,µ¯,ψ)∝ exp − − −
(σ2)n2 2 2σ2 2σ2
(cid:3) (cid:4) (cid:25)
τ 1
+2n 2 y¯tPtb ,
σ2
(cid:1) (cid:1)
where SST = ki=1 nj=i1 (yij − y¯)2 is the total sum of squares and y¯ =
(y¯ ,...,y¯ )t representsthesamplemeanvector.
1 k
UsingresultsinFerrándiz(1982)tointegrateoutψ andstandardintegration
methods for µ¯ and σ2, we obtain that the integrated likelihood results to be
thesameasthemarginallikelihood(2.2),showingthatnolossofinformation
derivesfromconsideringonlythemarginalexperiment,atleastwhenτ isthe
realparameterofinterest.
4 Choiceamongnoninformativepriors
Intheprevious section we have shown that no matter what reference prior is
used, the likelihood function for the parameter of interest τ turns out to be
alwaysthesame.Nowwediscusswhichreferencemarginalprio√rforτ should
be used, that is we want to compare the two priors π (τ) = 1/ τ(τ +2) and
d
492 F.Solarietal.
√
π (τ) = 1/ τ derivedfromthetwoapproachesdescribedintheprevioussec-
r
tions. Furthermore, a frequentist probability matching prior for τ can be also
derived;weillustratethisissueinSect.4.1.
4.1 Thefirst-ordermatchingpriorforτ
In this section we derive all the first-order matching priors for τ in the sense
that one-sided(cid:26)cred(cid:27)ible intervals and their frequentist coverage probabilities
agreeuptoO n−1 .FortheFisherinformationmatrixI of(τ,σ2,µ¯,ψ)given
in(3.1),itiseasytoshowthat
(cid:5)(cid:16) (cid:17) (cid:6)
I−1(τ,σ2,µ¯,ψ)= 2diag τ(τ +2) −τσ2 , σ2, 1(CCt)−1 .
n −τσ2 σ4 τ
Lete =(1,0,0,0(cid:7) )(cid:7);thene(cid:7)I−1e =2τ(τ +2)/nand
1 k−2 1 1
⎡√ ⎤
(cid:29) (cid:28)τ(τ +2)
I−1e 2 ⎢⎢− τ σ2⎥⎥
ν = (cid:28) 1 = ⎢ τ+2 ⎥.
e(cid:7)I−1e n⎣ 0 ⎦
1 1
0k−2
Following Dattaand Ghosh(1995),afirst-ordermatching priorπ(τ,σ2,µ¯,ψ)
forτ mustsatisfythedifferentialequation
(cid:29)
∂ (cid:30) ∂ τ
τ(τ +2)π(τ,σ2,µ¯,ψ)− σ2π(τ,σ2,µ¯,ψ)=0 (4.1)
∂τ ∂σ2 τ +2
√
Letπ∗ = τ(τ +2)σ2π(τ,σ2,µ¯,ψ).Clearly,(4.1)isthenequivalentto
∂ ∂
(τ +2) π∗−σ2 π∗ =0,
∂τ ∂σ2
which has a general solution π∗ = g((τ +2)σ2,µ¯,ψ), for some positive and
differentiablefunctiong(·,·,·).Thenageneralsolutionfor(4.1)is
1
π(τ,σ2,µ¯,ψ)= √ g((τ +2)σ2,µ¯,ψ).
τ(τ +2)σ2
It is easy to see that the reference prior π is a matching prior by choosing
(cid:31) d
g(·) = k−3 (sinψ)k−i−2, which is a constant on its first two arguments. On
i=1 i
theotherhand,π willnotsatisfy(4.1)andisnotamatchingprior.
r