Higher Accuracy for Bayesian and Frequentist Inference: Large Sample Theory for Small Sample Likelihood

0.35 MB·English

by  M. Bédard

StatisticalScience
2007,Vol.22,No.3,301–321
DOI:10.1214/07-STS240
(cid:13)c InstituteofMathematicalStatistics,2007
Higher Accuracy for Bayesian and
Frequentist Inference: Large Sample
Theory for Small Sample Likelihood
M. B´edard, D. A. S. Fraser and A. Wong
8
0
0 Abstract. Recent likelihood theory produces p-values that have remarkable accuracy
2 and wideapplicability. Thecalculations usefamiliar tools such as maximum likelihood

n values (MLEs), observed information and parameter rescaling. The usual evaluation
a of such p-values is by simulations, and such simulations do verify that the global dis-
J
tribution of the p-values is uniform(0, 1), to high accuracy in repeated sampling. The
4
derivation of the p-values, however, asserts a stronger statement, that they have a uni-
2
form(0, 1) distribution conditionally, given identified precision information provided

] by the data. We take a simple regression example that involves exact precision infor-
E
mation and use large sample techniques to extract highly accurate information as to
M
the statistical position of thedata pointwith respect to the parameter: specifically, we
.
t examine various p-values and Bayesian posterior survivor s-values for validity. With
a
t observed data we numerically evaluate the various p-values and s-values, and we also
s
[ record the related general formulas. We then assess the numerical values for accuracy

usingMarkov chain Monte Carlo(McMC) methods.We also proposesome third-order
1
likelihood-based procedures for obtaining means and variances of Bayesian posterior
v
1 distributions, again followed by McMC assessment. Finally we propose some adap-
5 tive McMC methods to improve the simulation acceptance rates. All these methods
7
are based on asymptotic analysis that derives from the effect of additional data. And
3
. the methods use simple calculations based on familiar maximizing values and related
1
informations.
0
8 The example illustrates the general formulas and the ease of calculations, while the
0 McMC assessments demonstrate the numerical validity of the p-values as percentage
:
v position of a data point. The example, however, is very simple and transparent, and
Xi thus gives little indication that in a wide generality of models the formulas do accu-
rately separate information for almost any parameter of interest, and then do give
r
a accurate p-value determinations from that information. As illustration an enigmatic
problem in the literature is discussed and simulations are recorded; various examples
in the literature are cited.
Key words and phrases: Asymptotics, Bayesian posterior s-value, canonical parame-
ter,defaultprior,higherorder,likelihood,maximumlikelihooddeparture,Metropolis–
Hastings algorithm, p-value, regression example, third order.
Myl`ene B´edard is Assistant Professor, D´epartement de
Math´ematiques et de Statistique, Universit´e de
York University, 4700 Keele Street, Toronto, Ontario,
Montr´eal, C.P. 6128, succ. Centre-ville, Montr´eal,
Canada M3J 1P3 e-mail: [email protected].
Qu´ebec, Canada H3C 3J7 e-mail:
[email protected]. D. A. S. Fraser is Professor, This is an electronic reprint of the original article
Department of Statistics, University of Toronto, 100 St published by the Institute of Mathematical Statistics in
George St., 6th floor, Toronto, Ontario, Canada M5S Statistical Science, 2007, Vol. 22, No. 3, 301–321. This
3G3 e-mail: [email protected]. A. Wong is reprint differs from the original in pagination and
Professor, Department of Mathematics and Statistics, typographic detail.
1
2 M. BE´DARD,D. A. S.FRASERAND A.WONG
1. INTRODUCTION Section11presentssomeintuitivethoughtsonthe
Metropolis–Hastings step in Markov chain Monte
We explore various large sample and likelihood
Carlo and then proposes several asymptotic and
methods for obtaining Bayesian and frequentist p-
adaptivemodificationsofthedirectMcMCapproach;
values from a regular statistical model and data.
theseareexploredforp-valuesinSections12and13.
Numerical values are obtained for a simple example
A controversial example is examined in Section 14,
that indicates the ease with which the methods can
and some concluding remarks are recorded in a dis-
be applied, given the typically available maximum
cussion Section 15.
likelihood and related calculations. The general for-
The Bayesian and frequentist methods give about
mulas are presented and discussed.
thesameanswerfortheexample.Infact,forthepar-
The example, ignoring the nonnormality of er-
ticular example they give theoretically the same an-
ror, is simple and transparent: one could plot the
swer,aconsequenceofthejudiciouschoiceofdefault
data, calculate means and standard deviations, do
prior for the Bayesian analysis. We do not address
the bootstrap, or even record likelihood and get here the manner of making such judicious choices
broadly about the same answer. But the large sam- or how the choice typically needs to be targeted on
ple techniques, more exactly data-accretion tech- the particular parameter of interest; this will be ad-
niques, provide accurate separation of component dressed subsequently.
parameter information, precisely summarize the
available information, and give accurate determina- 2. A SIMPLE EXAMPLE: DEPARTURE OF
tions of corresponding p-values. DATA FROM PARAMETER VALUE
InSection2wetakeapragmaticapproachandob-
Consider an example to illustrate the formulas
tain p-values using simple departure measures and
coming from large sample or, more exactly, data-
distributionalapproximationsrelated totheCentral
accretion techniques: a small data set involving a
LimitTheorem.TheninSection 3weformallyrefer-
response y with possible linear dependence on a re-
encethestatisticalmodelandobtainap-valuebased
lated variable x:
on the signed likelihood root. In Section 4 we add
x −3 −2 −1 0 1 2 3
a widely accepted default prior and obtain the pos- (1)
y −2.68 −4.02 −2.91 0.22 0.38 −0.28 0.03
terior survivor value, the analogue of the frequen-
tist p-value. For the example these require three- Theresponsevariability is taken to bethicker tailed
dimensional integration. than the normal, say the frequently suggested Stu-
But then in Section 5 we examine recently de- dent(7) distribution. Then with linear dependence
veloped likelihood-based approximations that have andconstantresponsevariability,wehavethemodel
third-order accuracy; numerical values are obtained 7
for the example, and the general formulas are dis- f(y;θ)=σ−7 h σ−1(y α βx ) ,
i i
{ − − }
cussed. In Section 6 we discuss the corresponding i=1
Y
frequentist third-order p-values. Numerical values where h(z) = Γ(4)/Γ(1/2)Γ(7/2)√7 (1 + z2/7)−4
for the example are then presented together with { }
is the Student(7) density. Or using quantile func-
various intermediate values that indicate how the
tion form, we can write y =α+βx +σz , where
i i i
calculations proceed.
thelatenterrorsz areindependentStudent(7).Now
i
In Sections 7 and 8 we consider exact p-values for suppose we are interested in assessing the response
the preceding methods, as derived by Markov chain dependence on x as given by the regression param-
Monte Carlo. As part of this we note that a N = eter β, with particular interest in whether β=1.
4 106 simulation in the particular context gives As background we note that the response data
×
about two-figure accuracy for probabilities, about were in fact generated from the given model with
thesameasthethird-orderapproximationmethods. Student(7) error and then rounded to two decimal
InSection 9webrieflydiscusstheroleof precision places; the parameter values used to generate the
information in the Bayesian and frequentist con- data were α= 0, β = 1 and σ =1. In passing we
texts. Section 10 looks directly at Bayesian means note that σ is an error scaling and does not record
and variances and how they can be approximated directly the error standard deviation. Also there is
by recent likelihood-based methods. Again Markov no implied connection between the number of ob-
chain Monte Carlo is used to evaluate the accuracy. servations n=7 and the degrees of freedom for the
HIGHERACCURACY FOR BAYESIANAND FREQUENTISTINFERENCE 3
error density df =7. The example has simple and p-value. Using the standard normal then as an ap-
transparentstructure,and weuseit to examine var- proximation, we obtain the approximate p-value
ious frequentist and Bayesian assessment methods; p =Φ(t0)=Φ( 1.445634)
N
we then apply Markov chain Monte Carlo sampling (3) −
to check the distributional validity of the resulting =0.07414=7.414%,
values.
which records the observed level of significance in
The example is simple and transparent and we
an elemental form, as just the percentage position
could probably do as well by plotting or by least
of the data point or the probability left of the data
squares and standard deviations. The likelihood point, under the hypothesis.
model theory, however, gives a precise separation of A simple modification hopefully to accommodate
information concerning parameters of interest and the estimation of error scaling is obtained by using
anaccurate determination of thep-valueor percent- the Student(5) distribution function, say H (), as a
5
·
age position of data with respect to a value for the revised approximation method. We then obtain the
parameter of interest. As a more general illustra- approximate p-value
tion, we later record simulation data for a challeng- p =H (t0)=H ( 1.445634)
S 5 5
ingexamplefromtheliterature.Wealsocitevarious (4) −
examples that have beenexamined in the literature. =0.10395=10.395%.
For theexampleconsidered,apragmaticfirststep
An alternative to the direct use of the large sample
is to use least squares to separate out the general
distribution theory is provided by the bootstrap ap-
location and the linear dependence on the related
proach. Using the least squares values, we separate
variable x. The fitted values for α and β are
the data values into a location or fit component yˆ
i
and a residual component y yˆ:
a= 1.322857, b=0.675000, i− i
− x −3 −2 −1 0 1 2 3
with residual length s = (SSR)1/2 = 2.660046 ob- yˆ −3.3479−2.6729−1.9979−1.3229−0.6479 0.0271 0.7021
tained from the sum of squares of residuals. y−yˆ 0.6679 −1.3471−0.9121 1.5429 1.0279 −0.3071−0.6721
The parameter β records the indicated linear de- For thebootstrapwe randomlysampletheresiduals
pendence of y on x. To assess the value β = 1 in with equal probability and add them back to the lo-
the presence of the data, we can reasonably exam- cation component, thus obtaining a bootstrap data
ine the raw departure b β= 0.325000, and then set,fromwhichwecalculatethebootstrapt-statistic
− −
standardizeitbyitsestimatedstandarddeviationto value
∗
obtain a standardized measure of departure of data ∗ b 1
t = − ,
from parameter value, giving s∗/√5√28
∗ ∗
b(y) β where b and s are the regression coefficient and
t= − ,
residual length from the bootstrap sample. We re-
s/√5√28
(2) peated this bootstrap sampling a convenient total
t0= 0.675000−1 = 1.445634. of N = 10,000 times, and the empirical distribu-
2.660046/√5√28 − tion function was evaluated at the observed t0 =
1.445634.Thisgaveanobservedbootstrapp-value:
To interpret this statistically we need information −
(5) p =Fˆ(t0)=0.1051=10.51%,
concerning the distribution of possible values for t BS
in the context where the true value of β is 1. To where Fˆ(t) is the empirical distribution function of
this end, some reference to large sample theory sug- the bootstrap t∗ values.
gests the use of the standard normal distribution With the present small sample size n=7 we can
function, say Φ(). The observed value of this dis- calculate the bootstrap p-value exactly, p , by
· ExBS
tribution function then gives an approximation to using equal probability for each of the 77=823,543
the percentage of possible values of t that would possible bootstrap samples. We then take
be less than the observed t0, in other words, to the p = proportion(t∗<t0)
ExBS
percentage position of the data with respect to the (6) {
hypothesized value β=1; this is called the observed +proportion(t∗ t0) /2,
≤ }
4 M. BE´DARD,D. A. S.FRASERAND A.WONG
Table 1 (9) =a 7logσ
Simple frequentist p-values for the regression example in −
Section 2 7 (y α βx )2
i i
4 log 1+ − − ;
− 7σ2
Measure of Distributional p-value i=1 (cid:26) (cid:27)
X
departure approximation β=1 β=1.5 β=2
theobservedlikelihoodL0(α,β,σ)andobservedlog-
t-statistic Normal 0.07414 0.03121 0.08189 likelihood ℓ0(α,β,σ) are obtained by substituting
t-statistic Student(5) 0.10395 0.02722 0.02100 the data y0 =(y0,...,y0) from the data array (1)
t-statistic Bootstrap N=104 0.10750 0.02790 0.03800 1 7
into the above expressions.
t-statistic Bootstrap (exact) N=770.10332 0.02833 0.03888
SLR Normal 0.05774 0.02148 0.04830 In a general context the observed likelihood func-
tion is given as
Ap-valuerecordsthepercentagepositionofthedatarelative
to a possible true value β for the parameter; β=1 is in fact L0(θ)=L(θ;y0)=cf(y0;θ)=f0(θ),
thetruevalueunderlyingthedata set; β=1.5 and β=2 are
other values that might have been of interest. Multiple zeros which is the observed density function, that is, the
are indicated by a superscript, thus0.02148=0.00148. statistical model f(y;θ) examined at the observed
data point y0; it records the amount of probability
called a mid-p-value, and obtain sitting at that data point, viewed as a function of
possible values for the parameter. The constant c is
(7) p =0.1033231=10.332%;
ExBS taken as arbitrary but positive and indicates that
with our particular rounded t0 value there were no only relative values from one θ value to another are
values at the boundary point. These four approx- of relevance given the data point. If we consider in
imate p-values make use of a pragmatic choice of generalhowlikelihooddependsondatawecanwrite
departure measure combined with distributional in-
L(θ;y)=cf(y;θ),
formation provided by Central Limit Theorem-type
approximations or by resampling from the nonpara- ℓ(θ;y)=a+logf(y;θ),
metric maximum likelihood distribution; and they
where c>0 and a,c are otherwise arbitrary for each
provide us with four determinations of the observed
choice of data point y.
percentagepositionofthedatarelativetothemodel
Suppose now that we are interested in a scalar
with β =1. Other initial departure measures could
componentψ=ψ(θ).Mostlikelihoodmethodsmake
have beenconsidered, as well as other distributional
use of maximum likelihood values (MLEs); we do,
approximations or determinations.
however, avoid referring to them as estimates, as
These p-values for testing β =1 are recorded in
they are useful but typically not directly as esti-
Table1togetherwithalikelihood-basedmethoddis-
mates. We write θˆ=argsupL(θ) for the value that
cussed in the next section. The table also records
maximizes L(θ). Also if ψ(θ) is a component
p-values for testing the more extreme β values, 1.5
and 2, with corresponding observed values t0 = parameter of particular interest, we write θˆψ =
3.669685 and t0= 5.893737. argsupψ(θ)=ψL(θ)forthevaluethatmaximizesL(θ)
− − subject to the interest parameter ψ(θ) having some
value ψ(θ)= ψ of special interest. Then based on
3. THE EXAMPLE: SIMPLE LIKELIHOOD
DEPARTURE MEASURE the likelihood L(θ) alone, an important departure
measure of data from ψ(θ)=ψ is obtained as the
More intrinsic departure measures are available
signed likelihood root (SLR)
from longavailable likelihood theory. Thelikelihood
functionin thepresentStudentregression context is (10) r =sign(ψˆ ψ)[2 ℓ(θˆ) ℓ(θˆ ) ]1/2;
ψ ψ
− { − }
L(α,β,σ;y) this measure examines how probability at the data
(8) point under the full model exceeds that when ψ(θ)
=cσ−7 7 1+ (yi−α−βxi)2 −4, is restricted to the value ψ, and theory has shown
7σ2 it to be fundamental. One way of viewing this mea-
i=1(cid:26) (cid:27)
Y sure is to picture how much the log-likelihood rises
with log-likelihood
from the maximum when ψ(θ)=ψ up to the over-
ℓ(α,β,σ;y) all maximum when θ is unrestricted, that is, ℓ(θˆ)
−
HIGHERACCURACY FOR BAYESIANAND FREQUENTISTINFERENCE 5
Fig. 1. For the example in Section 2: the third-order Bayesian survivor function, the third-order frequentist p-value and the
first-order SLR p-value.
ℓ(θˆ ), often written ℓˆ ℓ˜. We take this rise as hav- venient but can be awkward when the error distri-
ψ
−
ing quadratic form r2/2 in terms of some quantity butionitself hasdependenceontheparameter;forif
ψ
rψ; solve for rψ; and then attach the appropriate the distribution for the error itself has dependence
sign. For testing a true value ψ(θ)=ψ the signed ontheparameterratherthanasherebeingjustStu-
likelihood root rψ has to first order a standard nor- dent(7), then the maximum likelihood value could
mal limiting distribution, a follow-through from the
also have that parameter dependence and thus not
Central Limit Theorem. Related measures could be
be a statistic. From the preceding numerical values
basedon theslopeof thelog-likelihood atthetested
we obtain from (10) the signed likelihood root
value or on the displacement from θˆ to θˆ, but nei-
ψ
ther has the same mathematical invariance or the
(11) r = 1.574053.
β=1
same track record in applications. −
For many likelihood calculations, particularly re-
The corresponding observed p-value based on the
cent higher-order calculations, the computationally
first-order normal approximation for r is then
challenging aspects often arise in the maximization
steps rather than in other steps. (12) p =Φ(r )=5.774%;
SLR β=1
For our simpleregression example and the related
likelihood calculations we now use maximum like- this is also recorded in Table 1. In Figure 1 we
lihood rather than least squares variables; the ob- plot the SLR p-value against a full range of pos-
served values obtained by computer iteration are
sible β values; this is called the p-value function;
αˆ= 1.3504512, some related determinations discussed in later sec-
−
tionsarealsorecordedinthefigure.Otherlikelihood
βˆ=0.6504019,
departure measures based on the score and maxi-
σˆ=0.9641110 mumlikelihoodestimates aresometimes considered,
but they frequently have serious distributional and
for the overall MLEs, and
measurementbiasdifficulties.Bycontrast, theSLR-
α˜=αˆ = 1.366699,
β=1 based approximate p-value uses a departure mea-
−
β˜=βˆ =1, sure that directly relates to the statistical model;
β=1
it summarizes background information contained in
σ˜=σˆ =1.154527
β=1 themodelcombined with distributionalinformation
for the constrained MLEs when β =1. In passing derived from the large sample behavior of the like-
we note that the use of MLE variables can be con- lihood function.
6 M. BE´DARD,D. A. S.FRASERAND A.WONG
4. ANALYSIS WITH DEFAULT PRIOR this derives from parameter transformations on the
sample and parameter spaces and is referred to as
Analternativefirst-orderlikelihood-basedmethod
the left invariant prior (e.g., Jeffreys, 1946; Fraser,
comesfromtheuseofamodel-basedflatpriorcalled
1979):undertransformationsthatmakelocationand
adefaultpriortogether withconditional-probability
scale changes on the initial sample space, the dif-
type calculations, typically referred to as Bayes al-
ferential rewritten as, say, dα dβ dσ/σ3 remains
gorithm. Let π(θ) be a weight function for θ based
unchanged, is invariant. This left prior avoids the
on symmetries, invariance, or other relevant proper-
marginalization issues for certain parameter com-
ties of the model. The corresponding posterior dis-
ponents that are linear in a location parameteriza-
tribution viewed as providing inference information
tion impliedbyasymptotictheory (Fraser andReid,
concerning θ is
2002);formanyfamiliarparametersofinterest,how-
(13) π(θ y0)=cπ(θ)L(θ), ever, it can lead to the marginalization issues; and
|
furthermoreitdoesnotcorrespondtotheconfidence
where c now indicates the norming constant. This
theory pivotal inversion based on the usual equa-
default approach was implicit in Bayes (1763), was
tions (y α βx )/σ=z .
i i i
strongly promoted by Laplace (1812) and Jeffreys − −
For the parameter β, various approaches (e.g.,
(1946) and acquired the name inverse probability.
Dawid, Stone and Zidek, 1973; Fraser, 1979) sug-
More recently its development has been stimulated
gest the targeted prior
bytheValenciaconferences(see,e.g.,BayesianStatis-
tics 7, 2003); the default priors are called objective dα dβ dσ
π (θ)dθ= =dα dβ dlogσ,
β
priors but in such contexts the objective can only σ
refer to model structure, not to any objective fre-
and for many parameters having a certain linearity
quencies in the context being examined.
itdoes avoid themarginalization issue;somediscus-
For inference concerning an interest parameter,
sionofthislinearityandarelatedcurvaturemeasure
say ψ(θ), one might then reasonably calculate the
may be found in Fraser and Reid (2002); the curva-
marginal posterior density function
ture issue does not arise in the present problem for
theniceparametersα,β,σ andwedonotpursuethe
π(ψ y0)= π(θ y0)dλ,
| | issue here. From the transformation viewpoint this
Z
iscalledtherightinvariantprior.Thecorresponding
whereλisacomplementingnuisanceparameterhere
model-based posterior for θ is then given by
chosen so that θ is one–one equivalent to (ψ,λ) and
to have, say, Jacobian ∂θ/∂ψ∂λ 1 so that sup- π(α,β,τ;y0)dθ
| |≡
port volume corrections can be ignored.
(14) =cσ−7
This marginalization to obtain an inference dis-
tribution for a component parameter can produce 7 (y0 α βx )2 −4
misleading results (e.g., Dawid, Stone and Zidek, 1+ i − − i dαdβdlogσ;
· 7σ2
1973). To overcome this issue, a preferred approach iY=1(cid:26) (cid:27)
is to have a prior that depends on the particular forthis,ifwenowtakeθtobe(α,β,τ)=(α,β,logσ),
parameter of interest, and thus to use a targeted the implied prior is π(θ)=1 and it conforms to con-
prior πψ(θ)wherethesubscriptindicates thepartic- fidence inversion.
ular parameter being targeted (e.g., Jeffreys, 1946; In order then to obtain the marginal density for
Bernardo, 1979; Fraser et al., 2003). We do not ad- the interest parameter β, an integration over α and
dress this important issue of choosing default priors τ is required. Repeated numerical integration over
but do acknowledge that it is of major interest for twodimensionscanbequitefeasible butoften isnot
the Bayesian community at the present time and in easily implemented; we next consider some alterna-
part for the frequentist community. tive integration procedures.
For our simple regression example we might pos- For more general use of this Bayesian approach
sibly consider the model-based default prior π(θ) to the choice of the default prior becomes a crucial is-
be the invariant prior sue and is the focus of much present activity in the
Bayesian community; the term objective Bayesian
dα dβ dσ dα dβ dlogσ
π(θ)dθ= = ; prior is sometimes used in place of the term default
σ3 σ2
HIGHERACCURACY FOR BAYESIANAND FREQUENTISTINFERENCE 7
Bayesian prior, but this is misleading as the objec- The preceding marginal density can also be writ-
tive would indicate that it is describing the physical ten
context rather than being based as here on model ek/n L (ψ)
characteristics, a level removed from the physical π(ψ;y0)= P
(2π)d/2 L (ψˆ)
context. P
(18)
ˆj 1/2π(θˆ )
θθ ψ
5. THIRD ORDER WITH DEFAULT PRIOR | | ,
· j (θˆ ) π(θˆ)
(cid:26) λλ ψ (cid:27)
| |
Formanyregulardensities,posteriororotherwise,
where
the Laplace integration method provides an accu-
rate alternative route for obtaining the marginal L (ψ)=supL(ψ,λ)=L(ψ,λˆ )
P ψ
density of a component, such as β in our case. As λ
before, but in general notation, we use π(ψ,λ) for is the profile likelihood function for ψ, obtained by
the proposed prior and L(ψ,λ) for the likelihood. maximizing the full likelihood over λ for fixed value
Then with third-order accuracy, the marginal pos- ψ of the interest parameter ψ(θ).
terior density for ψ can be obtained by Laplace in- The methods inherent in the Laplace integration
tegration over the nuisance parameter divided by procedure can be described fairly easily. A regular
Laplace integration over the full parameter, giving function, here L(ψ,λ) for fixed ψ, whose logarithm
has additive and maximum likelihood value proper-
ek/n
π(ψ;y0)= ties under increasing sample size n can be rewritten
(2π)d/2
(15) f(λ)=eh(λ)
·e−rψ2/2 j |ˆj(θθθˆ| ) 1/2ππ((θˆθˆψ)), =cexp{−ˆjλλλ2/2}exp{aλ3/6n1/2+bλ4/24n}
(cid:26) λλ ψ (cid:27)
| | to third order as a function of λ, with obvious gen-
where k indicates a first-order constant: r2 is the eralization for vector λ; for this we are letting λ
ψ
likelihood ratio quantity designate a standardized departure of the original
(16) r2 =2 ℓ(θˆ) ℓ(θˆ ) λ from the maximizing value for f(λ) with ψ fixed.
ψ { − ψ } After expanding the second exponential in a power
discussed earlier but now used more generally with series and similarly for the log-prior, and then inte-
ψofdimensiond,nuisanceparameterλofdimension grating term by term with respect to λ, we obtain
m, and p=d+m; the Hessian matrices to fourth order,
∂2
f(λ)π(ψ,λ)dλ
j (θ;y)= ℓ(θ;y),
θθ −∂θ∂θ′ Z
∂2 =ek/n(2π)m/2 ˆjλλ(ψ)−1/2f(λˆψ)π(ψ,λˆψ),
j (θ;y)= ℓ(θ;y) | |
λλ −∂λ∂λ′ where ˆj is the negative Hessian with respect to λ as
are information functions for the full and for the evaluated at the maximum for the fixed ψ. The in-
nuisanceparameters,andhave dimensionsp pand tegrations are based on simple reference to the mul-
×
m m;they are justnegative second-derivative ma- tivariate normal integral; for some background, see
×
trices of the log-likelihood function, and when eval- Strawderman(2000)andforsomediscussionofterm
uated at θˆ and θˆ give the observed information by term integration, see Andrews, Fraser and Wong
ψ
matrices (2005).
For a scalar interest parameter ψ we can reason-
ˆj =j (θˆ;y),
θθ θθ ably be more interested in an integral of its den-
(17)
˜j =j (θˆ )=j (θˆ ;y). sity function, and particularly in the right tail in-
λλ λλ ψ λλ ψ
tegral called the posterior survivor function. Why
Numerically, the information matrices can typically the right tail? Consider the simple case of a vari-
be computed by taking second differences based on able x measuring a parameter ψ with error density
very small and equally spaced values for each coor- f(e) and distribution function F(e); we have: the
dinate. observed p-value is p0(ψ)=F0(ψ)=F(x0 ψ); the
−
8 M. BE´DARD,D. A. S.FRASERAND A.WONG
right tail posterior survivor function with a natu- Foroursimpleexampleandtestingβ=1,wehave
ral flat prior is s(ψ)= ∞f(x0;α)dα= ∞f(x0 r= 1.574053 from (11) and we have
ψ ψ − −
α)dα = F(x0 ψ); and these are equal. In more
R R
general model−situations the p-value as discussed in (22) qB= 0.9483686,
−
the next section records in a statistical sense where
where
the data value is with respect to ψ in a left-to-right
distributional sense and then corresponds as in the
ℓ (θˆ )= 5.868699,
simple case to the survivor posterior value s(ψ). β β=1 −
For the general case a highly accurate approxima-
and the full and constrained information matrices
tion to the posterior survivor function is available
for θ=(α,β,τ) are
from likelihood theory (see, e.g., Fraser, Reid and
Wu, 1999, generalizing DiCiccio and Martin, 1991):
ˆj =j (θˆ)
θθ θθ
s(ψ)=1 G(ψ y0) (23)
(19) − | 5.7894195 1.311060 −0.3286837
r =911 1.311060 27.288552 −1.3395559 ,
=Φ r r−1log =Φ(r∗). −0.3286837 −1.339556 12.1900132 !
− q
(cid:18) (cid:19)
Forthis,Gdesignatestheposteriordistributionfunc- (24) j (θˆ )=911 4.0240689 −0.3319537 ,
tion for ψ, r is the signed likelihood root r given λλ β −0.3319537 8.5925687
ψ (cid:18) (cid:19)
as (10) in Section 3, q is a score-type departure
B
with corresponding determinants
measure for ψ,
(20) qB=ℓψ(θˆψ) |ˆjλλˆj(θˆψ)| 1/2ππ((θˆθˆ)), (25) |ˆjθθ|=1892.702,
(cid:26) | θθ| (cid:27) ψ jλλ(θˆβ) =34.4669.
| |
where
∗
This gives r = 1.252169; the Bayesian posterior
∂ ∂
(21) ℓψ(θ)= ∂ψℓ(θˆψ)= ∂ψℓ(θ;y0) θˆψ survivor value fro−m (18) is then
(cid:12)
∗(cid:12)
is a score departure measure, and r (cid:12)is implicitly (26) sB(1)=0.1052542;
defined. Note that for convenience we have taken
′ ′
the full parameter θ to be given as (ψ,λ) in terms thisisrecordedinTable2togetherwiththeBayesian
of the components, and this applies as well to the survivor values for testing β=1.5 and β=2, as well
information matrices at (17). as some McMC validation results discussed in later
Also we have chosen to record the upper tail for sections.
presenting posterior probability from the posterior The first-order likelihood method in Section 3 re-
distribution; the interest parameter will often have
quires the full and the constrained maximum like-
physical meaning in a particular application and in- lihood values θˆ and θˆ with of course correspond-
ψ
vestigators will think in terms of a variable mea-
ing values for the log-likelihood function. In order
suring the parameter as, say, with a maximum like-
to take advantage of the approximate integration
lihood estimate. In such a framework the usual p-
formulas in this section, we require in addition the
valueis lefttail for thevariable andcorrespondingly
second-derivative values at each MLE; such deriva-
isrighttailfortheparameterasinthelocation case;
tives are of course also needed for familiar score and
accordingly for harmony between the two inference
MLE departure measures and typically can be ob-
approaches wetake thereference Bayesian probabil-
tained by differencing.
ity to be the upper tail survivor value. In the next
We can also calculate the Bayesian survivor value
section we will find that formula (20) using r from
(10) and q from (20) can be derived directly from s(β) for a range of values for β. For our special ex-
B
frequentist formulas in the next section by acting ample the Bayesian survivor function sB(β) is plot-
as if π(ψ;y0) in (15) were obtained from a location ted in Figure 1 together with the likelihood ratio
model π(ψ x;y0) with a nominal variable x taking p-value Φ(r ) and a third-order frequentist p-value
β
−
the observed value x=0. to be discussed in the next section.
HIGHERACCURACY FOR BAYESIANAND FREQUENTISTINFERENCE 9
6. THIRD-ORDER p-VALUE the density form
Recent likelihood methods give highly accurate (28) f¯(y;ϕ)dy= ek/n e−r2/2 ˆj 1/2dϕˆ
approximationsforfrequentistinferenceinmuchthe (2π)p/2 | ϕϕ|
same manner as for the Bayesian context just de-
ek/n
scribed. The methods require full and constrained (29) = e−r2/2 ˆj −1/2ds,
maximum likelihood values θˆ and θˆ , as well as (2π)p/2 | ϕϕ|
ψ
fullandconstrainedinformationdeterminants.They wherer2 isthelikelihoodratioquantityforassessing
also, however, need something more in the way of the full parameter ϕ and ˆj = (∂2/∂ϕ2)ℓ(ϕ;y0)
ϕϕ ϕˆ
− |
information from the model and data. The nature istheobservedinformationmatrixfromadatapoint
of this extra information can best be described in y. The renormalized (28) is third-order accurate.
terms of parameterization scaling or reexpression. For the importantspecial case of a scalar parame-
In particular, we need to express the initial parame- terϕ,acorrespondingdistributionfunctionapproxi-
terθ asacanonical-typeparameter,sayϕ=ϕ(θ).In mation was developed by Lugannani and
theBayesian context, such additional information is Rice(1980)andinanalternativeformbyBarndorff-
closely related to thedevelopment of an appropriate Nielsen (1991). Both versions use the signed likeli-
hood ratio r=r(ϕ,s) corresponding to (10), plus a
default prior; thus the use of a weighted likelihood
maximum likelihood value departure q=q(ϕ,s):
L(θ)π(θ) can partly be interpreted in terms of seek-
ingalocationreparameterization β=β(θ)suchthat q =(ϕˆ ϕ)ˆj1/2.
π(θ)dθ=dβ. Bayesian parameter reweightings have f − ϕϕ
The approximate distribution function in the
longbeensoughtinthedevelopmentalsequencefrom
Barndorff-Nielsen (1991) form for ϕˆ or s is then
invariant (Bayes, 1763; Laplace, 1812) to Jeffreys
(1946) to reference priors (Bernardo, 1979). (30) F¯(s;ϕ)=Φ r r−1log r =Φ(r∗),
In the frequentist context the accessible reparam- − q
(cid:18) f(cid:19)
eterizations areofanexponentialratherthanaloca-
andhasthird-orderaccuracy;thishasthesameform
tiontype,buttheygivesomeaccesstotherelatedlo-
as (19) but uses a different departure q appropriate
cation information. The exponential-type reparam-
to the present frequentist context. The similarity of
eterization can be examined initially in the context
the Bayesian formula (19) and the above frequen-
of an exponential model. To this effect, consider an
tist formula can appear more plausible by defining
exponential model with canonical parameter ϕ, the following reexpressions of the variable and the
′ parameter:
(27) f(y;ϕ)=exp ϕs(y) κ(ϕ) h(y),
{ − } ϕ
β(ϕ)= ˆj1/2dϕˆ,
where ϕ and s have the same dimension, say p. ϕϕ
Z
For such a model, the saddlepoint approximation s
b(s)= ˆj−1/2ds.
(Daniels,1954)canberemarkablyaccurate,andhas ϕϕ
Z
Table 2
For the simple regression example inSection 2, Bayesian s-values for assessing the values β=1,
β=1.5 and β=2: using the third-order formula (19); using the McMC (Section 8 with Normal
proposal); using AMcMC (Section 12 with adaptive choice of Student proposal)
Testprocedure s-value
β=1 β=1.5 β=2
Bayesian: third order (N=4×106) 0.10525 0.02725 0.03923
Bayesian: McMC (N=4×106, Normal) 0.10744 0.02841 0.02118
(Simulation SD) (0.03484) (0.03186) (0.04789)
{Acceptance rate} {41.9%} {41.9%} {41.9%}
Bayesian: AMcMC (N=4.106, adaptiveStudent) 0.10752 0.02836 0.02118
(Simulation SD) (0.03332) (0.03100) (0.04366)
{Acceptance rate} {51.1%} {51.1%} {51.1%}
10 M. BE´DARD,D. A. S.FRASERAND A.WONG
Interms of these,Welch and Peers (1963)showed in is the secondary departure measure. For this, θ =
effectthatbhasalocation modelf(b β)tosecond- (λ,ψ)hasbeenpresentedasacombinationofψ with
−
order accuracy. This has profound implications for a nuisance parameter λ which complements the in-
possible second order agreement between Bayesian terest parameter ψ; in the vector case we should
andfrequentist methodologies, butWelch and Peers perhaps write the combination in term of row vec-
′ ′ ′
(1963) presented their results in terms of the fre- tors as θ =(ψ ,λ). The scalar parameter χ(θ) is a
quentist approach of obtaining confidence bounds rotated coordinate of ϕ(θ) that acts as a surrogate
byintegrating likelihood withrespecttotheJeffreys for ψ(θ) and has linearity in terms of ϕ(θ). Explicit
(1946) prior i1/2(ϕ)dϕ. The same result with some formulas for the surrogate parameter χ(θ) and the
ϕϕ
greater generality can be obtained by Taylor series nuisance information are recorded in the Appendix;
expansion of an asymptotic model with scalar vari- somediscussionalsoappearsinthenextsection;the
ableandparameter; simplereexpressions ofvariable parentheses around (λλ) are to indicate that the in-
and parameter show that to second order the model formation has been recalibrated in terms of the new
canbewritteneitherasalocation orasanexponen- parameterization ϕ. All of this is easily accessible
tial model; see, for example, Cakmak et al. (1998) numerically, and uses primarily just the typical in-
andAndrews,Fraser andWong(2005).Asomewhat gredients of the Bayesian-type approximation.
similar result is available with vector variable and
parameter; see Cakmak, Fraser and Reid (1994). 7. THIRD ORDER FOR THE EXAMPLE
For the case of a location model f(s β) with a
− We have just described how third-order p-values
flat prior π(β)=c, the expression (20) for q sim-
B are available to assess a scalar parameter ψ(θ) in an
plifies to
exponential model. While exponential models are of
q =ℓ (β)ˆj−1/2. coursequiteimportant,theydorepresentaveryspe-
B β ββ
cialized type of model. However, recent likelihood
As ℓβ(β;s) = (∂/∂β)ℓ(β;s) = (∂/∂s)ℓ(s β) = theory has shown that for a general continuous sta-
− −
ℓ;s(β;s) from the location property, and ℓβ(β)= tistical model together with data, there exists a cor-
− −
ℓ;s(β,s)=ϕ from the exponential form (27), we ob- responding exponential model that provides highly
tain through simple algebra that qf =qB, which im- accurate third-order p-values for the original model
plies that the frequentist distribution function is and data, using the formulas in the preceding sec-
equal to the Bayesian survivor function, as would tion.
be expected. In other words, this Welch and Peers For our example in Section 3, the corresponding
(1963)Bayesian–frequentistequalityisobtainedfrom exponential model with data has the same observed
the Jeffreys prior for the scalar exponential model, log-likelihood ℓ(θ)=logf(y0;θ)given as(9)andhas
and thus as demonstrated by Cakmak et al. (1994, anominalreparameterization ϕ′=(ϕ ,ϕ ,ϕ )given
1 2 3
1998) has an extension for general asymptotic mod- as a row vector; some discussion is given later. See
els. The advantages of this scalar parameter use of also Davison, Fraser and Reid (2006). Thereparam-
the Jeffreys prior has recently been discussed for eterization is
thediscretebinomialdistributioncontextbyBrown,
′
ϕ(α,β,σ)
Cai and DasGupta (2001).
(32)
Now consider the vector exponential model (27)
7 (α+βx y0)/7σ2
with p-dimensional canonical parameter ϕ and p- =8 i− i (1,x ,d0),
1+(y0 α βx )2/7σ2 i i
dimensional canonical variable s, and suppose that Xi=1 i − − i
we are interested in a scalar component parameter
and is explained later in detail; here d0 is just the
ψ(ϕ) having reasonable smoothness properties.
standardized residual vector (y0 yˆ0)/σˆ0 recorded
Thesignedlikelihoodrootr=rψ in(10)isthepri- numerically preceding (35). The c−orresponding gen-
mary departure measure and a maximum likelihood
eral formulas arerecorded at the end of this section.
departure
Toobtainthep-valueforassessinganyscalarcom-
q (ψ)=sign(ψˆ ψ) ponent parameter, it suffices to treat the observed
f
−
(31) likelihood as a function of ϕ, which of course means
ˆjϕϕ 1/2 explicitly or implicitly that the observed informa-
χˆ χˆ | |
·| − ψ| j (θˆ ) tionsneedstobereexpressedorrecalibratedinterms
(cid:26) (λλ) ψ (cid:27)
| |