An Introduction to Bayesian Analysis: Theory and Methods - PDF Free Download (2024)

~ X 2 (1 + o(1))

by standard approximation to a normal tail (vide Feller (1973, p. 175) or Bahadur (1971, p. 1)). Thus~ logp --t -0 2 /2 and by (6.26), logB 01 --t-oo. This result is true as long as IX I > c(log n j n) 1 / 2 , c > V2. Such deviations are called moderate deviations, vide Rubin and Sethuraman (1965). Of course, even for such P-values, p rv (B 0 1/n) so that P-values are smaller by an order of magnitude. The conflict in measuring evidence remains but the decisions are the same. Ghosh et al. (2005) also pursue the comparison of the three measures of evidence based on the likelihood ratio, the P-value based on the likelihood ratio test, and the Bayes factor B 01 under general regularity conditions. 6.4.2 Pitman Alternative and Rescaled Priors

We consider once again the problem of testing H 0 : B = 0 versus H 1 : B -f. 0 on the basis of a random sample from N(B, 1). Suppose that the Pitman alternatives are the most important ones and the prior g 1 (B) under H 1 puts most of the mass on Pitman alternatives. One such prior is N(0,6jn). Then B 01 =

VS+i exp [- ~

(6

!

1) X

2

J.

If the P-value is close to zero, ..fiiiXI is large and therefore, B 01 is also close to zero, i.e., for these priors there is no paradox. The two measures are of the same order but the result of Berger and Sellke (1987) for symmetric unimodal priors still implies that P-value is smaller than the Bayes factor.

6.5 Bayesian P-value3 Even though valid Bayesian quantities such as Bayes factor and posterior probability of hypotheses are in principle the correct tools to measure the evidence for or against hypotheses, they are quite often, and especially in many practical situations, very difficult to compute. This is because either the alternatives are only very vaguely specified, vide (6.6), or very complicated. Also, in some cases one may not wish to compare two or more models but check how a model fits the data. Bayesian P-values have been proposed to deal with such problems. 3

Section 6.5 may be omitted at first reading.

180

6 Hypothesis Testing and Model Selection

Let M 0 be a target model, and departure from this model be of interest. If, under this model, X has density f(xlry), 17 E £, then for a Bayesian with prior 1r on 1], m1r(x) = J£ f(xlry)n(ry) dry, the prior predictive distribution is the actual predictive distribution of X. Therefore, if a model departure statistic

T(X) is available, then one can define the prior predictive P-value (or tail area under the predictive distribution) as

P = pm7C (T(X) 2': T(xobs)IMo), where Xobs is the observed value of X (see Box (1980)). Although it is true that this is a valid Bayesian quantity for model checking and it is useful in situations such as the ones described in Exercise 13 or Exercise 14, it does face the criticism that it may be influenced to a large extent by the prior 1r as can be seen in the following example.

Example 6.18. Let X 1 , X 2 , · · · , X n be a random sample from N ((), a 2 ) with both () and a 2 unknown. It is of interest to detect discrepancy in the mean of the model with the target model being M 0 : () = 0. Note that T = y'riX (actually its absolute value) is the natural model departure statistic for checking this. (a) Case 1. It is felt a priori that a 2 is known, or equivalently, we choose the prior on a 2 , which puts all its mass at some known constant a5. Then under M 0 , there are no unknown parameters and hence the prior predictive P-value is simply 2(1-

ex

(

~xT

) -n/2

,

which is an improper density, thus completely disallowing computation of the prior predictive P-value. (c) Case 3. Consider an inverse Gamma prior IG(v, j3) with the following density for a 2 : n(a 2 1v, j3) = (a 2 )-(v+l) exp(- ~) for a 2 > 0, where v and

A:)

j3 are specified positive constants. Because Tla 2 the predictive density ofT is then,

rv

N(O, a 2 ), under this prior

6.5 Bayesian P-value

181

Table 6.6. Normal Example: Prior Predictive P-values .5 .5 1 1 1 2 2 2 5 5 5 .5 1 .5 .5 1 1 2 1 2 2 2 p .300 .398 .506 .109 .189 .300 .017 .050 .122 .0001 .001 .011 v

.5

/3 .5

If 2v is an integer, under this predictive distribution,

Thus we obtain, P = pm~

~

(lXI

=

pm~

=

2(1-

lxobsiiMo)

(I_I_I > vnlxobsiiM) 0

Vffv - Vffv

p. ( vnlxobsl 2v

JNv

))

'

where Fzv is the c.d.f. of tzv· For foxobs = 1.96 and various values of v and {3, the corresponding values of the prior predictive P-values are displayed in Table 6.6. Further, note that p -t 1 as {3 -t oo for any fixed v > 0. Thus it is clear that the prior predictive P-value, in this example, depends crucially on the values of v and {3. What can be readily seen in this example is that if the prior 1r used is a poor choice, even an excellent model can come under suspicion upon employing the prior predictive P-value. Further, as indicated above, non-informative priors that are improper (thus making m7f improper too) will not allow computing such a tail area, a further undesirable feature. To rectify these problems, Gutman (1967), Rubin (1984), Meng (1994), and Gelman et al. (1996) suggest modifications by replacing 1r in m7f by n(rylxobs):

m*(xlxobs) = p* =

1

f(xlry)n(rylxobs) dry,

pm*([xobs)

(T(X)

~

and

T(Xobs)).

This is called the posterior predictive P-value. This removes some of the difficulties cited above. However, this version of Bayesian P-value has come

182

6 Hypothesis Testing and Model Selection

under severe criticism also. Bayarri and Berger (1998a) note that these modified quantities are not really Bayesian. To see this, they observe that there is "double use" of data in the above modifications: first to convert (a possibly improper) 1r(17) into a proper 7r(1JIXabs), and then again in computing the tail area of T(X) corresponding with T(xabs)· FUrthermore, for large sample sizes, the posterior distribution of 17 will essentially concentrate at ij, the MLE of ry, so that m*(xlxabs) will essentially equal f(xiiJ), a non-Bayesian object. In other words, the criticism is that, for large sample sizes the posterior predictive P-value will not achieve anything more than rediscovering the classical approach. Let us consider Example 6.18 again.

Example 6.19. (Example 6.18, continued.) Let us consider the non-informative prior 1r(a 2) IX 1/a2 again. Then, as before, because Tja 2 "'N(O, a 2), and

IX

n (-2 exp ( - 2a2 xobs

+ sobs 2 ))( a 2)-"'+1 2

'

Therefore, we see that, under the posterior predictive distribution,

Thus we obtain the posterior predictive P-value to be

where Fn is the distribution function of tn. This definition of a Bayesian P-value doesn't seem satisfactory. Let lxabsl --+ oo. Note that then p --+ 2 (1- Fn(vn)). Actually, p decreases to this lower bound as lxabsl--+ oo.

6.5 Bayesian P-value

Table 6.7. Values of Pn

183

= 2 (1- Fn(fo))

n 1 2 3 4 5 6 7 8 9 10 Pn .500 .293 .182 .116 .076 .050 .033 .022 .015 .010

Values of this lower bound for different n are shown in Table 6. 7. Note that these values have no serious relationship with the observations and hence cannot be really used for model checking. Bayarri and Berger (1998a) attribute this behavior to the 'double' use of the data, namely, the use of x in computing both the posterior distribution and the tail area probability of the posterior predictive distribution. In an effort to combine the desirable features of the prior predictive P-value and the posterior predictive P-value and eliminate the undesirable features, Bayarri and Berger (see Bayarri and Berger (1998a)) introduce the conditional predictive P-value. This quantity is based on the prior predictive distribution mrr but is more heavily influenced by the model than the prior. Further, noninformative priors can be used, and there is no double use of the data. The steps are as follows: An appropriate statistic U(X), not involving the model departure statistic T(X), is identified, the conditional predictive distribution m(tlu) is derived, and the conditional predictive P-value is defined as

Pc where Uabs (1998a).

= pm(.!uobs)

(T(X) 2': T(Xabs))'

U(xabs)· The following example is from Bayarri and Berger

Example 6.20. (Example 6.18, continued.) T = foX is the model departure statistic for checking discrepancy of the mean in the normal model. Let U(X) = s 2 = ~ 2:::7= 1 (X;- X) 2 . Note that nUia 2 "'a 2 x~_ 1 . Consider n(a 2 ) ex 1/a 2 again. Then n(a 2 IU = s 2 ) ex (a 2 )(n-l)/ 2 +1 exp( -ns 2 /(2a 2 )) is the density of inverse Gamma, and hence the conditional predictive density of T given s~bs is

1 t2 ) -n/2 ex ( 1+--2n 8 abs

Thus, under the conditional predictive distribution,

184

6 Hypothesis Testing and Model Selection

T Sobs

- - rv

tn-1,

and hence we obtain the conditional predictive P-value to be Pc

= pm(.ls~bs) (lXI 2': lxobsiMo) = pm(.ls~bsl (Jn=IIXI 2': Jn=IIxl IMo) Sobs

=

2 ( 1- Fn_ 1

(

Sobs

vn:~Xobsl)).

We have thus found a Bayesian interpretation for the classical P-value from the usual t-test. It is worth noting that s~bs was used to produce the posterior distribution to eliminate a 2 , and that Xobs was then used to compute the tail area probability. It is also to be noted that in this example, it was easy to find U(X), which eliminates a 2 upon conditioning, and that the conditional predictive distribution is a standard one. In general, however, even though this procedure seems satisfactory from a Bayesian point of view, there are problems related to identifying suitable U(X) and also computing tail areas from (quite often intractable) m(tluobs)· Another possibility is the partial posterior predictive P-value (see Bayarri and Berger (1998a) again) defined as follows: p*

=

pm*(.)

(T(X) 2': T(xobs)),

where the predictive density m* is obtained using a partial posterior density 7r* that does not use the information contained in tabs= T(xobs) and is given by m*(t)

=

1

fr(tlrJ)7r*(rJ) drJ,

with the partial posterior 1r* defined as 1r*(77) ex fxlr(Xobsltobs, rJ)7r(rJ)

fx(XobslrJ) ( ) ex fr(tobslrJ) 7r

7] •

Consider Example 6.18 again with 1r(a 2 ) ex 1/a2 . Note that because Xi are i.i.d. N(O, a 2 ) under M 0 ,

fx (xobs la 2 ) ex (a 2 ) -n/ 2 exp(- : 2 (x~bs 2

ex so that

+ s~bs))

fx (Xobs la 2 ) ( a 2 ) -(n-I)/ 2 exp(-

: 2

2

s~bs),

6.6 Robust Bayesian Outlier Detection

f XIX_(Xobs ~-Xobs, a 2)

185

2 ). ex (a 2)-(n-l)/2 exp ( - -n-2 sobs 2a

Therefore, 1f* (a2)

ex (a2) -(n-l)/2+1 exp(- 2:2 s~bs),

and is the same as 1r(a 2 ls~bJ, which was used to obtain the conditional predictive density earlier. Thus, in this example, the partial predictive P-value is the same as the conditional predictive P-value. Because this alternative version p* does not require the choice of the statistic U, it appears this method may be used for any suitable goodness-of-fit test statistic T. However, we have not seen such work.

6.6 Robust Bayesian Outlier Detection Because a Bayes factor is a weighted likelihood ratio, it can also be used for checking whether an observation should be considered an outlier with respect to a certain target model relative to an alternative model. One such approach is as follows. Recall the model selection set-up as given in (6.1). X having density f(xiB) is observed, and it is of interest to compare two models Mo and M 1 given by Mo : X has density j(xiB) where 8 E 8o; M1 :X has density j(xiB) where 8 E 8 1 .

For i = 1, 2, gi(B) is the prior density of B, conditional on Mi being the true model. To compare Mo and M 1 on the basis of a random sample x = (x 1 , ... , Xn) the Bayes factor is given by

mo(x) B01 (x) = -(-), m1x where mi(x) = fe, f(xiB)gi(B) dB fori = 1, 2. To measure the effect on the Bayes factor of observation xd one could use the quantity kd = log

(

B(x) )

B(x-d)

,

(6.27)

where B(x-d) is the Bayes factor excluding observation xd. If kd < 0, then when observation xd is deleted there is an increase of evidence for Mo. Consequently, observation Xd itself favors model M 1 . The extent to which xd favors M 1 determines whether it can be considered an outlier under model M 0 . Similarly, a positive value for kd implies that xd favors M 0 . Pettit and Young (1990) discuss how kd can be effectively used to detect outliers when the prior is non-informative. The same analysis can be done with informative priors also. This assumes that the conditional prior densities g 0 and g 1 can be fully

186

6 Hypothesis Testing and Model Selection

specified. However, we take the robust Bayesian point of view that only certain broad features of these densities, such as symmetry and unimodality, can be specified, and hence we can only state that g0 or g 1 belongs to a certain class of densities as determined by the features specified. Because kd, derived from the Bayes factor, is the Bayesian quantity of inferential interest here, upper and lower bounds on kd over classes of prior densities are required. We shall illustrate this approach with a precise null hypothesis. Then we have the problem of comparing

Mo : e = eo versus M 1 : e -1- eo using a random sample from a population with density f(xle). Under M 1 , suppose e has the prior density g, g E T. The Bayes factors with all the observations and without the dth observation, respectively, are

B 9(x)

= foo~oa f(xle)g(e)de'

f(x-dleo) B9(x-d) = fo#a f(x-dle)g(e) de. Because f(xle) = f(xdle)j(x-dle), we get kd 9 =log

'

f(xleo) fo#a f(x-dle)g(e) del . [ f(x-dleo) .ftOofOa f(xle)g(e) de

l

foo~oa f(xle)g( e) de =log f(xdleo)- log [ fo#a f(x-dle)g(e) de .

(6.28)

Now note that to find the extreme values of kd, 9 , it is enough to find the extreme values of (6.29)

over the set

r. Further, this optimization problem can be rewritten as follows: sup hd 9 = sup

9EG

,

foo~oo f(xdle)J(±-dle)g(e) de

~--'-';;-------::-c,--------,--::.,----;--::-:---:c:--

fo#(Jo f(±-dle)g(e) de

9EG

sup 9*EG*

.

mf hd 9

9EG

r

loo~oa

f(xdle)g*(e) de,

(6.30)

foo~oa f(xdle)J(±-dle)g(e) de

.

= mf ---'---';;---,---------,----,---:--:---

'

9EG

=

inf 9*EG*

fo#Oo f(±-dle)g(e) de

r

loo~oa

f(xdle)g*(e) de,

(6.31)

6.6 Robust Bayesian Outlier Detection

187

where

G

*

{ =

*

g(B)f(;r_diB)

*( )

g :g B

=

}

fuo~oag(u)f(;r_dlu)du'gEG .

Equations (6.30) and (6.31) indicate how optimization of ratio of integrals can be transformed to optimization of integrals themselves. Consider the case where r is the class A of all prior densities on the set {B : B -=f. B0 }. Then we have the following result. Theorem 6.21. If f(;r-diB) > 0 for each B -=f. Bo, sup hd,g = sup f(xdiB), and

gEA

(6.32)

Oo/=Oo

in£ hd,g = in£ f(xdiB).

gEA

(6.33)

Ooj=Oo

Proof. From (6.30) and (6.31) above, sup hd, 9 =

gEA

sup

g* EA *

{

f(xdiB)g*(B) dB,

jl}oj=l} 0

where

A* =

{ * *( ) g :g B

g(B)f(;r_diB)

=

}

fuo~oog(u)f(;r_dlu)du'gEA .

Now note that extreme points of A* are point masses. Proof for the infimum is similar. The corresponding extreme values of kd are sup kd,g =log f(xdiBo) -log in£ f(xdiB),

(6.34)

in£ kd,g = logf(xdiBo) -log sup f(xdiB).

(6.35)

gEA

Ooj=Oo

gEA

l}oj=l} 0

Example 6.22. Suppose we have a sample of size n from N(B, a 2 ) with known a 2 • Then, from (6.34) and (6.35), sup kd,g = - 1 2 sup [ (xd- B) 2 - (xd- Bo) 2] gEA 2a Oo/=Oo = oo, and in£ kd,g =

gEA

~ in£ 2a O#Oo

[(xd- B) 2

-

(xd- Bo) 2 ]

(xd- Bo) 2 2a 2 It can be readily seen from the above bounds on kd that no observation, however large in magnitude, will be considered an outlier here. This is because A is too large a class of prior densities.

188

6 Hypothesis Testing and Model Selection

Instead, consider the class G of all N(B 0 , 7 2 ) priors with 7 2 > 76 > 0. Note that 7 2 close to 0 will make M 1 indistinguishable from M 0 , and hence it is necessary to consider 7 2 bounded away from 0. Then for g E G

where g* is the density of N(m, 82 ) with

Note, therefore, that hd,g = hd,g(xd) is just the density of N(m, a 2 evaluated at Xd· Thus,

+ 82 )

For each xd, one just needs to graphically examine the extremes of the expression above as a function of 7 2 to determine if that particular observation should be considered an outlier. Delampady (1999) discusses these results and also results for some larger nonparametric classes of prior densities.

6. 7 Nonsubjective Bayes Factors 4 Consider two models M 0 and M 1 for data X with density fi(xiOi) under model Mi, Oi being an unknown parameter of dimension Pi, i = 0, 1. Given prior specifications gi(Oi) for parameter Oi, the Bayes factor of M 1 to M 0 is obtained as B = m1(x) = f h(xiOI)gl(0 1)d0 1 (6.36) 10

mo(x)

J fo(xl0o)9o(Oo)d0o ·

Here mi(x) is the marginal density of X under Mi, i = 0, 1. When subjective specification of prior distributions is not possible, which is frequently the case, one would look for automatic method that uses standard noninformative priors. 4

Section 6.7 may be omitted at first reading.

6. 7 Nonsubjective Bayes Factors

189

There are, however, difficulties with (6.36) for noninformative priors that are typically improper. If gi are improper, these are defined only up to arbitrary multiplicative constants ci; cigi has as much validity as gi. This implies that (cl/c0 )B 10 has as much validity as B 10 . Thus the Bayes factor is determined only up to an arbitrary multiplicative constant. This indeterminacy, noted by Jeffreys (1961), has been the main motivation of new objective methods. We shall confine mainly to the nested case where fo and h are of the same functional form and fo(xiOo) is the same as fr(xjOI) with some of the co-ordinates of 0 1 specified. However, the methods described below can also be used for non-nested models. It may be mentioned here that use of diffuse (flat) proper prior does not provide a good solution to the problem. Also, truncation of noninformative priors leads to a large penalty for the more complex model. An example follows.

Example 6.23. (Testing normal mean with known variance.) Suppose we observe X= (X 1 , ... ,Xn)· Under Mo,Xi are i.i.d. N(0,1) and under M 1 , Xi are i.i.d. N(B, 1), BE R is the unknown mean. With the uniform noninformative prior gf (B) c for B under M 1 , the Bayes factor of M 1 to M 0 is given by B~ = ~cn- 1 1 2 exp[nX 2 /2].

=

e ::;

If one uses a uniform prior over - K ::; K, then for large K, the new Bayes factor B{b is approximately 1/(2Kc) times Bfo. Thus for large K, one is heavily biased against M 1 . This is reminiscent of the phenomenon observed by Lindley (1957). A similar conclusion is obtained if one uses a diffuse proper prior such as a normal prior N(O, T2 ), with variance T2 large. The corresponding Bayes factor is

which is approximately (nT 2 )- 1 12 exp[nX 2 /2] for large values ofnT 2 and hence can be made arbitrarily small by choosing a large value of T 2 . Also Bp0orm is highly non-robust with respect to the choice of T2 , and this non-robustness plays the same role as indeterminacy. The expressions for Bfo and Bp0orm clearly indicate similar behavior of these Bayes factors and the similar roles of V2ifc and (T 2 + 1/n)- 1/2. A solution to the above problem with improper priors is to use part of the data as a training sample. The data are divided into two parts, X= (X 1, X 2 ). The first part X 1 is used as a training sample to obtain proper posterior distributions for the parameters (given X I) starting from the noninformative priors

190

6 Hypothesis Testing and Model Selection

These proper posteriors are then used as priors to compute the Bayes factor with the remainder of the data (X 2 ). This conditional Bayes factor, conditioned on X 1 , can be expressed as

(6.37) where mi(X 1) is the marginal density of X 1 under Mi, i = 0, 1. Note that if the priors cigi, i = 0, 1, are used to compute B 10 (X 1), the arbitrary constant multiplier cl/ c0 of B 10 is cancelled by (co/ cl) of m0 (X I) /m 1 (X 1) so that the indeterminacy of the Bayes factor is removed in (6.37). A part of the data, X 1 , may be used as a training sample as described above if the corresponding posteriors 9i(BiiX I), i = 0,1 are proper or, equivalently, the marginal densities mi(X I) of X 1 under Mi, i = 0,1 are finite. One would naturally use minimal amount of data as such a training sample leaving most part of the data for model comparison. As in Berger and Pericchi (1996a), a training sample X 1 may be called proper if 0 < mi(X 1 ) < oo, i = 0, 1 and minimal if it is proper and no subset of it is proper.

Example 6.24. (Testing normal mean with known variance.) Consider the setup of Example 6.23 and the uniform noninformative prior g 1 (B) = 1 for B under M 1 . The minimal training samples are subsamples of size 1 with mo(Xi) = (1/J27r)e-xU 2 and m1(Xi) = 1. Example 6.25. (Testing normal mean with variance unknown.) Let X (X1, ... ,Xn)· Mo : X1, ... , Xn are i.i.d. N(O, a5), M1: X1, ... ,Xn are i.i.d. N(JL,ai). Consider the noninformative priors go(ao) = 1/ao under Mo and 91(JL,al) = 1/a1 . Here m 1 (Xi) = oo for a single observation Xi and a minimal training sample consists of two distinct observations xi, xj and for such a training sample (Xi, Xj), (6.38)

6. 7.1 The Intrinsic Bayes Factor

As described above, a solution to the problem with improper priors is obtained using a conditional Bayes factor B 10 (X 1), conditioned on a training sample X 1 . However, this conditional Bayes factor depends on the choice of

6.7 Nonsubjective Bayes Factors

191

the training sample X 1 . Let X(l), l = 1,2, ... L be the list of all possible minimal training samples. Berger and Pericchi (1996a) suggest considering all these minimal training samples and taking average of the corresponding L conditional Bayes factors B 10 (X(l))'s to obtain what is called the intrinsic Bayes factor (IBF). For example, taking an arithmetic average leads to the arithmetric intrinsic Bayes factor (AIBF) L

AIBF 10

=

_!_"'"""'

B 10

mo(X(l)) L ~ ml(X(l))

(6.39)

and the geometric average gives the geometric intrinsic Bayes factor (GIBF)

GIBFw = Bw (

IT mm

1

(X(l))) (X(l))

112

,

(6.40)

the sum and product in (6.39) and (6.40) being taken over the L possible training samples X(l), l = 1, ... , L. Berger and Pericchi (1996a) also suggest using trimmed averages or the median (complete trimming) of the conditional Bayes factors when taking an average of all of them does not seem reasonable (e.g., when the conditional Bayes factors vary much). AIBF and GIBF have good properties but are affected by outliers. If the sample size is very small, using a part of the sample as a training sample may be impractical, and Berger and Pericchi (1996a) recommend using expected intrinsic Bayes factors that replace the averages in (6.39) and (6.40) by their expectations, evaluated at the MLE under the more complex model M 1 . For more details, see Berger and Pericchi (1996a). Situations in which the IBF reduces simply to the Bayes factor B 10 with respect to the noninformative priors are given in Berger et al. (1998). The AIBF is justified by the possibility of its correspondence with actual Bayes factors with respect to "reasonable" priors at least asymptotically. Berger and Pericchi (1996a, 2001) have argued that these priors, known as "intrinsic" priors, may be considered to be natural "default" priors for the testing problems. The intrinsic priors are discussed here in Subsection 6. 7.3.

6. 7.2 The Fractional Bayes Factor O'Hagan (1994, 1995) proposes a solution using a fractional part of the full likelihood in place of using parts of the sample as training samples and averaging over them. The resulting "partial" Bayes factor, called the fractional Bayes factor (FBF), is given by

FBFw = ml(X,b) mo(X, b) where b is a fraction and

192

6 Hypothesis Testing and Model Selection

Note that F BF10 can also be written as

where (6.41) To make FBF comparable with the IBF, one may take b = mjn where m is the size of a minimal training sample as defined above and n is the total sample size. O'Hagan also recommends other choices of b, e.g., b = ynjn or lognjn.

We now illustrate through a number of examples.

Example 6.26. (Testing normal mean with known variance.) Consider the setup of Example 6.23. The Bayes factor with the noninformative prior g1 (B) = 1 was obtained as Bw =

J'iiin- 1 12 exp[n.X 2 /2] = J'iiin- 1 12 Aw

where A10 is the likelihood ratio statistic. Bayes factor conditioned on Xi is

Thus n

AIBFw

= n- 1

LB

n

10 (Xi)

= n- 3 / 2 exp(nX 2 /2)

i=1

GIBFw

L exp( -Xl/2), i=1

=

n- 1 / 2 exp[n.X 2 /2- (1/2n)

L xf].

The median IBF (MIBF) is obtained as the median of the set of values Bw(Xi), i = 1, 2, · · ·, n. The FBF with a fraction o < b < 1 is

F BF10

=

b1 / 2 exp[n(1 - b)X 2 /2]

= n- 112 exp[(n- 1)X 2 /2],

if b = 1/n.

Example 6.27. (Testing normal mean with variance unknown.) Consider the setup of Example 6.25. For the standard noninformative priors considered in this example, we have

6.7 Nonsubjective Bayes Factors B 10 -

r(n-1)

X

2

X

TG) 1

AI BF1o = B10 x -(n) 1f

r(n-1)

FBF10 =

v;r/(-B')

[

x

~n

1

("""' X2)n/2 Ut

[~(Xi- X) 2]Y "'"'

L.....t

IXi- Xjl

2 1::;;i<jC::n

~: (x:

193

(X 2

+

"

x2 ]~-1 , -"X) 2

X 2) J

2 with b = -:;;,·

Example 6.28. (normal linear model.) This example is from Berger and Pericchi (1996b, 2001 ). Berger and Pericchi determined the IBF for linear models for both the nested and non-nested case. We consider here only the nested case. Suppose for the data Y(n x 1) we consider the linear models

where /3i = (f3i 1, Pi2, · · ·,Pip,)' and a} are unknown parameters, and Xi is an n x Pi known design matrix of rank Pi < n. Consider priors of the form

Here qi = 0 gives the reference prior of Berger and Bernardo (1992a), and qi = Pi corresponds with the Jeffreys prior. For the nested case, when Mo is nested in M 1, Berger and Pericchi ( 1996b) consider a modified Jeffreys prior for which qo = 0 and q1 = P1 -Po· The integrated likelihoods mi(Y) with these priors can be obtained as

where C is a constant not depending on i, and Ri is the residual sum of squares under Mi, i = 0, 1. The Bayes factor B10 with the modified Jeffreys prior is then given by B 10 = (27r)(Pl-Po)/2

I

X' X 11/2 0 0 1 2

IX~X1I /

(R ) _Q

(n-po)/2

R1

( 6.42)

Also, one can see that a minimal training sample Y(l) in this case is a sample of size m = p 1+ 1 such that for the corresponding design matrices Xi ( l) (under Mi), x:(l)Xi(l), i = 0,1 are nonsingular. The ratio m 0(Y(l))/m 1(Y(l)) can be obtained from the expression of B 10 by inverting it and replacing n, X 0 , X1, Ro, and R1 by m, Xo(l), X 1(l), Ro(l), and R1(l), respectively, where Ri (l) is the residual sum of squares corresponding to (Y (l)) under Mi, i = 0, 1. Thus the conditional Bayes factor B 10 (Y(l)), conditioned on Y(l) is given by B10(Y(l))

IX~Xol1/2 (Ro)(n-po)/21X~(l)X1(l)l1/2 (R1(l))(Pl-Po+l)/2 IX~X111/2

R1

IX~(l)Xo(Z)I1/2

Ro(l)

194

6 Hypothesis Testing and Model Selection

One may now find an average of these conditional Bayes factors to find an IBF. For example, an arithmetic mean of B 10 (Y(l))'s for all possible minimal training samples Y(l)'s gives the AIBF, and a median gives the median IBF. In case of fractional Bayes factor, one obtains that (see, for example, Berger and Pericchi, 2001, page 152), with mf(X) as defined in (6.41), mg(X)

=

(_!!__) (Pl -po)/2

mt(X)

with b =

J

X~ X ljl/2 ( Rl) (m-po)/2

jX~Xojl/2

27r

Ro

m/n and hence FBF10 = b(Pl-Pol/2(Ro/R 1 )(n-m)/2.

See also O'Hagan (1995) in this context. For more examples, see Berger and Pericchi (1996a, 1996b, 2001) and O'Hagan (1995). Several other methods have been proposed as solutions to the problem with noninformative priors. Smith and Spiegelhalter (1980) and Spiegelhalter and Smith (1982) propose the imaginary minimal sample device; see also Ghosh and Samanta (2002b) for a generalization. Berger and Pericchi (2001) present comparison of four methods including the IBF and FBF with illustration through a number of examples. Ghosh and Samanta (2002b) discuss a unified derivation of some of the methods that shows that in some qualitative and conceptual sense, these methods are close to each other.

6. 7.3 Intrinsic Priors Given a default Bayes factor such as the IBF or FBF, a natural question is whether it corresponds with an actual Bayes factor based on some priors at least approximately. If such priors exist, they are called intrinsic priors. A default Bayes factor such as IBF can then be calculated as an actual Bayes factor using the intrinsic prior, and one need not consider all possible training samples and average over them. A "reasonable" intrinsic prior that corresponds to a naturally developed good default Bayes factor may be considered with be a natural default prior for the given testing or model selection problem. On the other hand, a particular default Bayesian method may be evaluated on the basis of the corresponding intrinsic prior depending on how "reasonable" the intrinsic prior is. Berger and Pericchi (1996a) describe how one can obtain intrinsic priors using an asymptotic argument. We begin with an example.

Example 6.29. (Example 6.26, continued.) Suppose that for some proper prior

1r(O) under model M 1 , BFf0

~

AI BF10

(6.43)

6.7 Nonsubjective Bayes Factors

195

where BF!0 denotes the Bayes factor based on a prior n(O) under M 1 . Using Laplace approximation (Section 4.3) to the integrated likelihood under M 1 , we have

BF1f 10

e

~ h(XIO)

fo(XIO = 0)

n- 1 12 n(e)v2,;i(deti)- 1 12

e

where is the MLE of under M 1 , and f is the observed Fisher information number. Thus using the expression for the AIBF in this example, and noting that i = 1, (6.43) can be expressed as 1

n(O) ~ (1jv2,;i)n

n

2:: exp( -Xl/2). i=l

As the RHS converges to (1/v'21f)E0 (e-xU 2 ) = (1/v'21f)(1/-/2)e-li probability one under any this suggests the intrinsic prior

e,

2 /

4

with

which is a N(O, 2) density. One can easily verify that

BF!0 / AI BFw -+ 1 with probability one under any e, i.e., the AIBF is approximately the same as the Bayes factor with an N(O, 2) prior fore under M 1 . If one considers the FBF, one can directly show that the FBF, with fraction b, is exactly equal to the Bayes factor with a N(O, (b- 1 - 1)/n) prior. Let us now consider the general case. Let B 10 be the Bayes factor of M 1 to = 0, 1. We illustrate below with the AIBF. Treatment for the other IBFs and FBF will be similar. Recall that

M 0 with noninformative priors gi(Oi) for ()i under Mi, i

L

1 '""'mo(X(l)) AIBFw = BwB01 where B01 = - L ( ( )) .

L

l=l

m1 X l

Suppose for some priors Jri under Mi, i = 0, 1, AIBF10 is approximately equal to the Bayes factor based on no and n 1 , denoted BF10 (n 0 , nl). Using Laplace approximation (Section 4.3) to both the numerator and denominator of B 10 (see 6.36), AIBFw can be shown to be approximately equal to (6.44) where n denotes the sample size, Pi is the dimension of ()i, Oi is the MLE of ()i, and ii is the observed Fisher information matrix under Mi, i = 0, 1. The same approximation applied to BF10 (n 0 , nl), yields the approximation

196

6 Hypothesis Testing and Model Selection

h (XIOdn1 (01) (2n jn )Pii 2li1l- 112 fo (XIOo)no (Oo) (2n jn )Pol 2 lio l- 1/ 2

(6.45)

to BF10 (n 0 , nl). We assume that conditions for the Laplace approximation hold for the given models. To find the intrinsic priors, we equate (6.44) with (6.45) and this yields (6.46) Berger and Pericchi (1996a) obtain the intrinsic prior determining equations by taking limits on both sides of (6.46) under M 0 and M 1. Assume that, as n-+ oo, under M1, 01-+ 61, Oo-+ a(Ol), and Bo1-+ Br(Ol); under Mo, Oo-+ Oo, 0 1 -+ b(Oo), and Bo1-+ B 0(0a). The equations obtained by Berger and Pericchi (1996a) are

When M 0 is nested in M 1, Berger and Pericchi suggested the solution (6.48) However, this may not be the unique solution to (6.47). See also Dmochowski (1994) in this context. Example 6.30. (Example 6.27, continued.) A solution to the intrinsic prior determining equations suggested by Berger and Pericchi (see (6.48)) is (6.49) where

"'xi

where Z1 = (X1-X2) 2j(2ar) and Z2 = (X1 +X2) 2j(2ai) "'noncentral 2 x with d.f. = 1 and noncentrality parameter >. = 2J1 2far. Also, Z 1 and Z 2 are independent. Using the representation of a noncentral x2 density as an (infinite) weighted sum of central x2 densities, we have

6.7 Nonsubjective Bayes Factors E [ z1 where

wj

rv

Zf ] + z2

=

~ e->-./2 (>../2)i E

L-.. ]=0

XI+2j and is independent of

[

j!

zl·

z1

Zf

+ wj

]

197 (6.50)

We then have

and the intrinsic priors are given by

with It is to be noted that f~oo 1r1 (f.LirYI)df.L = 1.

Example 6.31. (Testing normal mean with variance unknown.) This is from Berger and Pericchi (1996a). Consider the setup of Example 6.25 with the same prior g0 under M 0 but in place of the standard noninformative prior 9I(f.L,rYI) = 1/rYI use the Jeffreys prior gi(f.L,rYI) = 1/rYi. In this case, a minimal training sample consists of two distinct observations Xi, Xi for which

Proceeding as in the previous example, noting that mo(X1,X2) m1(X1,X2) where zl and z2 are as above, and using (6.50), the intrinsic priors are obtained as

Here n 1 (f.LjrY1) is a proper prior, very close to the Cauchy (0, rY 1) prior for f.L, which was suggested by Jeffreys (1961) as a default proper prior for f.L (given rY1); see Subsection 2.7.2.

Example 6.32. Consider a negative binomial experiment; Bernoulli trials, each having probability B of success, are independently performed until a total of n successes is accumulated. On the basis of the outcome of this experiment we want to test the null hypothesis H 0 : B = ~against the alternative H 1 : B # ~­ We consider this problem as choosing between the two models M 0 : B = ~ andM1=BE(0,1).

198

6 Hypothesis Testing and Model Selection

The data may be looked upon as n observations X 1 , ... , Xn where X 1 denotes number of failures before the first success, and for i = 2, · · ·, n, Xi denotes number of failures between (i - 1)th success and ith success. The random variables X 1 , .•• , Xn are i.i.d. with a common geometric distribution with probability mass function

P(Xi=x) =Bx(l-8), x=0,1,2, ... The likelihood function is

f(x 1, · · ·, X n le) -_ eE;'= 1 Xi (1 _ B)n · We consider the Jeffreys prior

which is improper. The Bayes factor with this prior is

Minimal training samples are of size 1, and the AIBF is given by

AIBFlO

= BlO

1 X :;:;

n

L

2Xi

+1

2Xi+2 .

i=l

Let

B*(B) = E [2X1 + 1] = ~ (2x + 1) ex( 1 _B). e 2 x, +2 L 2x+2 x=O

Then the intrinsic prior is 1

n(B) =

00

e- 1 12 (1- e)- 1 B*(B) = 4 L(2x + 1)Txex- 112 • x=O

Simplification yields

n(B)

=

(B- 112

+ 8 112 /2)(2- e)- 2 .

We now consider a simple example from Lindley and Phillips (1976), also presented in Carlin and Louis (1996, Chapter 1). In 12 independent tosses of a coin, one observes 9 heads and 3 tails, the last toss yielding a tail. It is shown that one gets different results according to a binomial or a negative binomial likelihood. Let us consider the problem of testing the null hypothesis H 0 : B = 1/2 against the alternative H 1 : B =f. 1/2 where B denotes the probability of head in a trial. If a binomial model is assumed, the random observable X is the number of heads observed in a fixed number of 12 tosses. One rejects H 0 for large values of the statistic IX- 6!, and the corresponding

6.8 Exercises

199

P-value is 0.150. On the other hand, if a negative binomial model is assumed, the random observable X is the number of heads before the third trial appears. Note that expected value of X under H 0 is 3. Suppose one rejects H 0 for large values of IX- 3J. Then the corresponding P-value is 0.0325. Thus with the usual 5% Type 1 error level, the two model assumptions lead to different decisions. Let us now use a Bayes test for this problem. For the binomial model, the Jeffreys prior is proportional to e- 112 (1 - e)- 112 , which can be normalized to get a proper prior. For the negative binomial model, the data can be treated as three i.i.d. geometrically distributed random variables, as described above. The Bayes factor under the binomial model (with Jeffreys prior) and the Bayes factor under the negative binomial model (with the intrinsic prior) are respectively 1.079 and 1.424. They are different as were the P-values of classical statistics, but unlike the P-values, the Bayes factors are quite close.

6.8 Exercises 1. Assume a sharp null and continuity of the null distribution of the test statistic. (a) Calculate EHo (P-value) and EH 0 (P-valueJP-value < a), where 0 < a < 1 is the Type 1 error probability. (b) In view of your answer to (a), do you think 2(P-value) is a better measure of evidence against H 0 than P-value? 2. Suppose X rv N(B, 1) and consider the two hypothesis testing problems: H 0 : B = -1 versus H 1 : B = 1; H~

: B = 1 versus H; : B = -1.

Find the Bayes factor of H 0 relative to H 1 and that of H 0 relative to Hi_ if (a) x = 0 is observed, and (b) x = 1 is observed. Compute the classical P-value in both cases. 3. Refer to Example 6.3. Take r = 20', but keep the other parameter values unchanged. Compute B 01 for the same values oft and n as used in Table 6.1. 4. Suppose X rv N(B, 1) and consider testing

H0 : B = 0 versus H 1 : B -=/= 0. For three different values of x, x = 0, 1, 2, compute the upper and lower bounds on Bayes factors when the prior on B under the alternative hypothesis lies in (a) rA ={all prior distributions on R},

(b)

TN=

{N(O,r 2 ),r 2 > 0},

(c) Fs ={all symmetric (about 0) prior distributions on R},

200

6 Hypothesis Testing and Model Selection

(d) rsu ={all unimodal priors on R, symmetric about 0}. Compute the classical P-value for each x value. What is the implication of rN c rsu c rs c rA? 5. Let X"' B(m, B), and let it be of interest to test H0

: ()

=

1

versus H1 :():f.

2

1

2·

If m = 10 and observed data is x = 8, compute the upper and lower bounds on Bayes factors when the prior on () under the alternative hypothesis lies in (a) rA ={all prior distributions on (0, 1)}, (b) rs ={Beta( a, a), a> 0}, (c) rs = {all symmetric (about ~) priors on (0, 1)}, (d) rsu = {all unimodal priors on (0, 1), symmetric about ~ }. Compute the classical P-value also. 6. Refer to Example 6.7. 1 (a) Show that B(GA,x) = exp(- ~ ), P(HolGA, x) = [1+ 1 -;;:o exp(~ . (b) Show that, if t:::; 1, B(Gus,x) = 1, and P(HolGus,x) =?To. (c) Show that, ift:::; 1, B(GNonx) = 1, and P(HoJGNonx) =?To. Ift > 1, B(GNor,x) = texp(-(t 2 -1)/2). 7. Suppose XJO has the tp(3, 0, Ip) distribution with density

)r

J(xJO) ex ( 1 +

8. 9. 10. 11.

1

3

) -(3+p)/2

(x- O)'(x- 0)

,

and it is of interest to test H 0 : 0 = 0 versus H 1 : 0 :f. 0. Show that this testing problem is invariant under the group of all orthogonal transformations. Refer to Example 6.13. Show that the testing problem mentioned there is invariant under the group of scale transformations. In Example 6.16, find the maximal invariants in the sample space and the parameter space. In Example 6.17, find the maximal invariants in the sample space and the parameter space. Let X I() rv N ((), 1) and consider testing Ho:

JB- Bol

::::: 0.1 versus H1 :

JB- Bol > 0.1.

Suppose x = () 0 + 1.97 is observed. (a) Compute the P-value. (b) Compute Bo 1 and P(Holx) under the two priors, N(B 0 , T 2 ), with (0.148) 2 and U(Bo- 1, 80 + 1). 12. Let XJp"' Binomial(10, p). Consider the two models:

M0 : p =

1

2

versus M1 : p :f.

1

2.

T

2

=

6.8 Exercises

201

Under M 1 , consider the following three priors for p: (i) U(O, 1), (ii) Beta(10, 10), and (iii) Beta(100, 10). If four observations, x = 0, 3, 5, 7, and 10 are available, compute kd given in Equation (6.27) for each observation, and for each of the priors and check which of the observations may be considered outliers under M 0 . 13. (Box (1980)) Let X 1 , X 2 , · · ·, Xn be a random sample from N (e, a- 2 ) with both and cr 2 unknown. It is of interest to detect discrepancy in the variance of the model with the target model being

e

Mo: a- 2 =

a-5,

and

e "'-' N(f-1, T 2 ),

where f-1 and T 2 are specified. (a) Show that the predictive distribution of (X1, X2, · · ·, Xn) under Mo is multivariate normal with covariance matrix a5In + T 2 11' and E(Xi) = f-1, for i = 1, 2, ... , n. (b) Show that under this predictive distribution,

(c) Derive and justify the prior predictive P-value based on the model departure statistic T(X). Apply this to data, x = (8, 5, 4, 7), and a-5 = 1, f-1 = 0, T 2 = 2. (c) What is the classical P-value for testing H 0 : a- 2 = a-5 in this problem? 14. (Box (1980)) Suppose that under the target model, for i = 1, 2, ... , n,

Yii,Bo, 0, a- 2 = ,Bo + x~O + Ei, Ei "'-' N(O, cr 2 ) i.i.d., ,6olcr 2 "'-' N(JLo,ccr 2 ),0ia 2 "'-' Nv(Oo,cr 2 F), a- 2 "'-'inverse Gamma( a,')'), where c, f-Lo, 0 0 , r, a and I' are specified. Assume the standard linear regression model notation ofy = ,60 1+X0+E, and suppose that X'l = 0. Further assume that, given a- 2 , conditionally ,60 , 0 and E are independent. Also, let /30 and 0 be the least squares estimates of ,60 and 0, respectively, and RSS = (y- /!Jol- XO)'(y- 0 1- XO). (a) Show that under the target model, conditionally on cr 2 , the predictive density of y is proportional to

/3

( a-2)-n/2

exp(--l-((/3o-JJ,o)2 +RSS 2cr2

c+l/n

(b) Prove that the predictive distribution of y under the target model is a multivariate t. (c) Show that the joint predictive density of ( RS S, 0) is proportional to

202

6 Hypothesis Testing and Model Selection {

21'

~

~

+ RSS + (8- 8o)'((X' x)- 1 + r- 1 ) - 1 (8- 8 0 )

}-(n+a-1)/2

(d) Derive the prior predictive distribution of

(0- 8o)'((X' x)- 1 + r- 1 ) - 1 (0- 8 0 ) T(y)

21'

=

+ RSS

.

(e) Using an appropriately scaled T(y) as the model departure statistic derive the prior predictive P-value. 15. Consider the same linear regression set-up as in Exercise 14, but let the target model now be

Mo: 8

=

O,,Bo[a 2

,...,

N(f.1o,ca 2 ),a 2

,...,

inverse Gamma(a:,--y).

Assuming/' to be close to 0, use ~I

T( ) y

=

~

8 X'X8 RSS

as the model departure statistic to derive the prior predictive P-value. Compare it with the classical P-value for testing H 0 : 8 = 0. 16. Consider the same problem as in Exercise 15, but let the target model be

Mo: 8

=

~I

~

2

2

2

1

O,,Bo[a ,..., N(f.lo,ca ),n(a) ex 2· a

Using T(y) = 8 X' X8 as the model departure statistic and RSS as the conditioning statistic, derive the conditional predictive P-value. Compute the partial predictive P-value using the same model departure statistic. Compare these with the classical P-value for testing H 0 : 8 = 0. 17. Let X 1 , X2, · · ·, Xn be i.i.d. with density

j(x[-\,0)

=

-\exp(--\(x- O)),x > 0,

where ,\ > 0 and -oo < 0 < oo are both unknown. Let the target model be

Mo : 0 = 0, n(-\) ex

1

>;·

Suppose the smallest order statistic, T = Xc 1l is considered a suitable model departure statistic for this problem. (a) Show that T[-\,..., exponential(n-\) under M 0 . (b) Show that -\[xabs ,..., Gamma(n, nxabs) under Mo. (c) Show that m (t IXabs ) = (

-n

nxabs

t

+ Xabs )n+ 1

(d) Compute the posterior predictive P-value. (e) Show that as tabs -----+ oo, the posterior predictive P-value does not necessarily approach 0. (Note that tabs :::; Xabs -----+ oo also.)

6.8 Exercises

203

18. (Contingency table) Casella and Berger (1990) present the following twoway table, which is the outcome of a famous medical experiment conducted by Joseph Lister. Lister performed 75 amputations with and without using carbolic acid. Patient Lived? Yes No

Carbolic Acid Used? Yes No 34 19 16 6

Test for association of patient mortality with the use of carbolic acid on the basis of the above data using (a) BIC and (b) the classical likelihood ratio test. Discuss the different probabilistic interpretations underlying the two tests. 19. On the basis of the data on food poisoning presented in Table 2.1, you have to test whether potato salad was the cause. (Do this separately for Crab-meat and No Crab-meat). (a) Formulate this as a problem of testing a sharp null against the alternative that the null is false. (b) Test the sharp null using BIC. (c) Test the same null using the classical likelihood ratio test. (d) Discuss whether the notions of classical Type 1 and Type 2 error probabilities make sense here. 20. Using the BIC analyze the data of Problem 19 to explore whether crabmeat also contributed to food poisoning. 21. (Goodness of fit test). Feller (1973) presents the following data on bombing of London during World War II. The entire area of South London is divided into 576 small regions of equal area and the number (nk) of regions with exactly k bomb hits are recorded.

Test the null hypothesis that bombing was at random rather than the general belief that special targets were being bombed. (Hint: Under H 0 use the Poisson model, under the alternative use the full multinomial model with 5 parameters and use BIC.) 22. (Hald's regression data). We present below a small set of data on heat evolved during the hardening of Portland cement and four variables that may be related to it (Woods et al. (1932), pp. 635-649). The sample size (n) is 13. The regressor variables (in percent of the weight) are x 1 = calcium aluminate (3Cao.Ab0 3), x 2 = tricalcium silicate (3CaO.Si0 2), x 3 = tetracalcium alumino ferrite (4Ca0.Ab0 3.Fe203), and X4 = dicalcium

204

6 Hypothesis Testing and Model Selection

Table 6.8. Cement Hardening Data X1 X2 X3 X4

7 1 11 11 7 11 3 1 2 21 1 11 10

26 6 29 15 56 8 31 8 52 6 55 9 7117 31 22 54 18 47 4 40 23 66 9 68 8

60 52 20 47 33 22 6 44 22 26 34 12 12

Y 78.6 74.3 104.3 87.6 95.9 109.2 102.7 72.5 93.1 115.9 83.8 113.3 109.4

silicate (2CaO.Si0 2 ). The response variable is y = total calories given off during hardening per gram of cement after 180 days. Usually such a data set is analyzed using normal linear regression model of the form

where pis the number of regressor variables in the model, {30 , {31 , ... {3p are unknown parameters, and Ei 's are independent errors having a N(O, a 2 ) distribution. There are a number of possible models depending on which regressor variables are kept in the model. Analyze the data and choose one from this set of possible models using (a) BIC, (b) AIBF of the full model relative to all possible models.

7 Bayesian Computations

Bayesian analysis requires computation of expectations and quantiles of probability distributions that arise as posterior distributions. Modes of the densities of such distributions are also sometimes used. The standard Bayes estimate is the posterior mean, which is also the Bayes rule under the squared error loss. Its accuracy is assessed using the posterior variance, which is again an expected value. Posterior median is sometimes utilized, and to provide Bayesian credible regions, quantiles of posterior distributions are needed. If conjugate priors are not used, as is mostly the case these days, posterior distributions will not be standard distributions and hence the required Bayesian quantities (i.e., posterior quantities of inferential interest) cannot be computed in closed form. Thus special techniques are needed for Bayesian computations.

Example 7.1. Suppose X is N(B,rr 2 ) with known rr 2 and a Cauchy(f1,T) prior one is considered appropriate from robustness considerations (see Chapter 3, Example 3.20). Then

and hence the posterior mean and variance are

Note that the above integrals cannot be computed in closed form, but various numerical integration techniques such as IMSL routines or Gaussian quadrature can be efficiently used to obtain very good approximations of these. On the other hand, the following example provides a more difficult problem.

206

7 Bayesian Computations

Example 7. 2. Suppose X 1 , X 2 , ... , X k are independent Poisson counts with Xi '""Poisson(Oi)· ei are a priori considered related, and a joint multivariate normal prior distribution on their logarithm is assumed. Specifically, let vi = log(Oi) be the ith element of v and suppose

where 1 is the k-vector with all elements being 1, and JL, constants. Then, because

and n(v) ex exp (-

2 ~ 2 (v -J-Ll)' ((1- p)Ik + pll')-

T

1

2

and pare known

(v -J-Ll))

we have that n(v[x) ex exp {- 2:::::7= 1 { ev; - vixi}-

~(v- J-Ll)' ((1- p)Ik + pll')- 1 (v -J-Ll)}.

Therefore, if the posterior mean of ej is of interest, we need to compute

) E 1f(n·[ u1 x

= exp {- 2:::::7= 1 { ev' -

= E1f(

( ·)[ ) exp v1 x

= fnk exp(vj)g(v[x) dv

JR.k g (V [X ) d V

where g(v[x)

Vixi}-

2

;2 (v- J-Ll)' ((1- p)h + pll')-

1

(v- J-Ll)}.

This is a ratio of two k-dimensional integrals, and as k grows, the integrals become less and less easy to work with. Numerical integration techniques fail to be an efficient technique in this case. This problem, known as the curse of dimensionality, is due to the fact that the size of the part of the space that is not relevant for the computation of the integral grows very fast with the dimension. Consequently, the error in approximation associated with this numerical method increases as the power of the dimension k, making the technique inefficient. In fact, numerical integration techniques are presently not preferred except for single and two-dimensional integrals. The recent popularity of Bayesian approach to statistical applications is mainly due to advances in statistical computing. These include the E-M algorithm discussed in Section 7.2 and the Markov chain Monte Carlo (MCMC) sampling techniques that are discussed in Section 7.4. As we see later, Bayesian analysis of real-life problems invariably involves difficult computations while MCMC techniques such as Gibbs sampling (Section 7.4.4) and MetropolisHastings algorithm (M-H) (Section 7.4.3) have rendered some of these very difficult computational tasks quite feasible.

7.1 Analytic Approximation

207

7.1 Analytic Approximation This is exactly what we saw in Section 4.3.2 where we derived analytic large sample approximations for certain integrals using the Laplace approximation. Specifically, suppose

E1f (g(O)[x)

=

Ink g(O)J(x[B)Ir(B) d() Ink f(x[O)Jr(O) d()

is the Bayesian quantity of interest where g, 8. First, consider any integral of the form

J,

and

(7.1)

are smooth functions of

1r

I= { q(O)exp(-nh(O)) dO, Jnk where h is a smooth function with -h having its unique maximum at 8. Then, as indicated in Section 4.3.1 for the univariate case, the Laplace method involves expanding q and h about Bin a Taylor series. Let h' and q' denote the vectors of partial derivatives of h and q, respectively, and Llh and Llq denote the Hessians of hand q. Then writing

h(O)

~

1

~

~

~

h(O)

=

h(O) + -(0- 0)' Llh(O)(O- 0) + · · · and

~

q(O)

~

1

~

~

~

2

~

=

~

+ (0- O)'h'(O) + 2(0- 0)' Llh(O)(O- 0) + · · ·

=

q(O)

~

~

1

~

~

~

+ (0- O)'q'(O) + 2(0- O)'Llq(O)(O- 0) + · · ·,

we obtain

I=

Lk {

q(O) + (0- O)'q'(O) +

~(()- 0)' Llq(O)(O- 0) +

... }

~ ~ xe-n h(B) exp ( - n (0- ~ 0)' Llh(O)(O0) + · · ·) d()

2

=

e-nh(0)( 27r)k/2n-ki2[Llh(O)[-l/2 { q(O)

+ O(n-1)},

(7.2)

which is exactly (4.16). Upon applying this to both the numerator and denominator of (7.1) separately (with q equal tog and 1), a first-order approximation

easily emerges. It also indicates that a second-order approximation may be available if further terms in the Taylor series expansion are retained in the approximation. Suppose that gin (7.1) is positive, and let -nh(O) = logf(x[O)+logJr(O), -nh*(O) = -nh(O) + logg(O). Now apply (7.2) to both the numerator and

208

7 Bayesian Computations

denominator of (7.1) with q equal to 1. Then, letting()* denote the maximum of -h*, E = Llh" 1 (0), E* = Llh"}(o;), as mentioned in Section 4.3.2, Tierney and Kadane (1986) obtain the fantastic approximation

(7.3)

which they call fully exponential. This technique can be used in Example 7.2. Note that to derive the approximation in (7.3), it is enough to have the probability distribution of g(O) concentrate away from the origin on the positive side. Therefore, often when g is non-positive, (7.3) can be applied after adding a large positive constant to g, and this constant is to be subtracted after obtaining the approximation. Some other analytic approximations are also available. Angers and Delampady (1997) use an exponential approximation for a probability distribution that concentrates near the origin. We will not be emphasizing any of these techniques here, including the many numerical integration methods mentioned previously, because the availability of powerful and all-purpose simulation methods have rendered them less powerful.

7.2 The E-M Algorithm We shall use a slightly different notation here. Suppose YIO has density f(yiO), and suppose the prior on () is 1r(O), resulting in the posterior density 1r(Oiy). When 1r(Oiy) is computationally difficult to handle, as is usually the case, there are some 'data augmentation' methods that can help. The idea is to augment the observed data y with missing or latent data z to obtain the 'complete' data x = (y, z) so that the augmented posterior density 1r(Oix) = 1r(Oiy,z) is computationally easy to handle. The E-M algorithm (see Dempster et al. (1977), Tanner (1991), or McLachlan and Krishnan (1997)) is the simplest among such data augmentation methods. In our context, the E-M algorithm is meant for computing the posterior mode. However, if data augmentation yields a computationally simple posterior distribution, there are more powerful computational tools available that can provide a lot more information on the posterior distribution as will be seen later in this chapter. The basic steps in the iterations of the E-M algorithm are the following. Let p(zly, 0) be the predictive density of Z given y and an estimate 0 of 0. ')

A

(

i)

A

(

i)

Find z' = E(Ziy, () ), where() is the estimate of() used at the ith step of the iteration. Note the similarity with estimating missing values. Use z(i) to (

augment y and maximize 1r(Oiy,z') to obtain() ')

(

A

(i+I)

. Then find z

('

t+I

)

using

and continue this iteration. This combination of expectation followed by maximization in each iteration gives its name to the E-M algorithm.

O(i+IJ

7.2 The E-M Algorithm

209

Implementation of the E-M Algorithm Note that because 1f(Oiy) = 1f(0, zly)/p(zly, 0), we have that log1f(Oiy) = log1f(0, zly) -logp(zly, 0). Taking expectation with respect to ZIO(i), y on both sides, we get log1f(Oiy) = =

J

J

A(i) log1f(0, zly)p(zly, 0 ) dz-

A(i) logp(zly, O)p(zly, 0 ) dz

Q(O, O(i))- H(O, O(i))

(7.4)

(where Q and H are according to the notation of Dempster et al. (1977)). Then, the general E-M algorithm involves the following two steps in the ith iteration: i) E-Step: Calculate Q(O, 0 ); A

(

A

M-Step: Maximize Q(O, 0

(

i)

A

(

) with respect to 0 and obtain 0

i)

such that

maxQ(O,O(i)) = Q(o,o(i)). (}

Note that log 7r(O(i+l) IY) -log 1f( o(i) IY) = { Q(O(i+l)' o(i)) - Q(O(i)' o(i))} - { H(O(i+ll, O(i))- H(O(i), O(i))}. From the E-M algorithm, we have that Q(O 2 + 1 , 0 for any 0, ,

)

A

(

(

H(O, O(i))- H(O(i), O(i))

J -J [ _J

J

=

A(i) logp(zly, O)p(zly, 0 ) dz-

-

p(zly, 0) ] ( I ll(i)) dz log p(zly, O(i)) p z y, u

--

log

i)

) ?: Q(O

logp(zly, 0

() 2

( 2')

A

,

(

i)

). Further,

A(i) )p(zly, 0 ) dz

[p(zly, O(i))] (i) p(zly, O) p(zly, 0 ) dz

::; 0, because, for any two densities P1 and P2, Jlog(pl(x)/P2(x))Pl(x)dx is the Kullback-Leibler distance between p 1 and p 2 , which is at least 0. Therefore,

H( O(i+l), O(i)) - H( O(i), O(i)) S 0, and hence 1f(0(i+l)IY)?: 1f(0(i)IY) for any iteration i. Therefore, starting from any point, the E-M algorithm can usually be expected to converge to a local maximum.

210

7 Bayesian Computations

Table 7.1. Genetic Linkage Data Cell Count Probability Y1 = 125 2+4

~: : ~~ 1~~ =:j

Y4

=

34

4

f!_

Example 7.3. (genetic linkage model.) Consider the data from Rao (1973) on a certain recombination rate in genetics (see Sorensen and Gianola (2002) for details). Here 197 counts are classified into 4 categories as shown in Table 7.1, along with the corresponding theoretical cell probabilities. The multinomial mass function in this example is given by J(y[O) ex: (2 + O)Y1 (1- O)Y 2 +Y3 0Y4 , so that under the uniform(O, 1) prior one, the observed posterior density is given by

This is not a standard density due to the presence of 2 +e. If we split the first cell into two with probabilities 1/2 and e14, respectively, the complete data will be given by x = (x1, x2, x3, x4, xs), where x1 + x2 = Yl, X3 = y2, X4 = Y3 and x 5 = Y4· The augmented posterior density will then be given by

which corresponds with the Beta density. The E-step of E-M consists of obtaining Q(O, {j(i))

=

E [ (X2 + X 5 ) loge+ (X3 + X4) log(1 - 0) [y, {j(i)]

= {

E

[x2[y, {j(i)]

+ Y4} log 0 + (y2 + Y3) log(1- 0).

(7.5)

TheM-step involves finding {j(i+l) to maximize (7.5). We can do this by solving JLQ(O {j(i)) = 0 so that

ae

'

'

{j( i+l) =

E

E

[x

[x2[y, {j(i)] + Y4

-~-"---~,--------=----

[x2[y, {j(i)]

[x

+ Y4 + Y2 + Y3

Now note that E 2 [y, {j(i)] = E 2 [X 1 + X 2 , {j(i)], and that '( ·l e(i) /4 X 2 [X 1 + X 2, 0' "'binomial(X1 + X2, l/HB(i)/ 4 ). Therefore,

and hence

(7.6)

7.3 Monte Carlo Sampling

211

Table 7.2. E-M Iterations for Genetic Linkage Data Example

Iteration i B(i) 1 .60825 2 .62432 3 .62648 4 .62678 5 .62682 6 .62682

t)(i)

'(i+l) _

e

-

+ Y4 + Y2 + Y3 + Y4

Y1 2+B<')

--~e<~,~)~-----------

Yl 2+B<')

(7.7)

In our example, (7. 7) converges to iJ = .62682 in 5 iterations starting from fj(o) = .5 as shown in Table 7.2.

7.3 Monte Carlo Sampling Consider an expectation that is not available in closed form. An alternative to numerical integration or analytic approximation to compute this is statistical sampling. This probabilistic technique is a familiar tool in statistical inference. To estimate a population mean or a population proportion, a natural approach is to gather a large sample from this population and to consider the corresponding sample mean or the sample proportion. The law of large numbers guarantees that the estimates so obtained will be good provided the sample is large enough. Specifically, let f be a probability density function (or a mass function) and suppose the quantity of interest is a finite expectation of the form

Ejh(X) =

i

(7.8)

h(x)f(x) dx

(or the corresponding sum in the discrete case). If i.i.d. observations xl' x2, ... can be generated from the density j, then

(7.9) converges in probability (or even almost surely) to Ejh(X). This justifies using hm as an approximation for Ejh(X) for large m. To provide a measure of accuracy or the extent of error in the approximation, we can again use a statistical technique and compute the standard error. If VarJh(X) is finite, then VarJ(hm) = VarJh(X)/m. Further, Var 1 h(X) can be estimated by

=

E 1 h 2 (X)-

(E1 h(X))

2

212

7 Bayesian Computations

and hence the standard error of

hm

can be estimated by

If one wishes, confidence intervals for Eth(X) can also be provided using the central limit theorem. Because

Vm (hm- Eth(X)) Sm

(

---+ N 0,1

)

m-HxJ

in distribution, (hm - Zo:j2Sm/ yrn, hm + Zo:j2Sm/ yrn) can be used as an approximate 100(1- o:)% confidence interval for Eth(X), with zo:; 2 denoting the 100(1 - o:/2)% quantile of standard normal. The above discussion suggests that if we want to approximate the posterior mean, we could try to generate i.i.d. observations from the posterior distribution and consider the mean of this sample. This is rarely useful because most often the posterior distribution will be a non-standard distribution which may not easily allow sampling from it. Note that there are other possibilities as seen below.

Example

7.4.

(Example 7.1 continued.) Recall that

1

2

E,.(Bix)

=

J~00 Bexp (- ((l;~) ) ( T2 + (B- J.L) 2) -- dB (J

f~oo exp (- (02~~)2) (T2 + (B- J.L)2)-1 dB 1

f~oo B { ~¢ (~)} ( T2 + (B- J.L) 2) -- dB f~oo { ~¢ (o;;x)} (T2 + (B- J.L)2)-1 dB , where¢ denotes the density of standard normal. Thus E,.(Bix) is the ratio of expectation of h(B) = Bj(T 2 + (B- J.L) 2) to that of h(B) = 1/(T 2 + (B- J.L) 2), both expectations being with respect to the N(x, CT 2 ) distribution. Therefore, we simply sample B1 ,B2 , ... from N(x,CT 2 ) and use ~ E,.(Bix)

"'m

B· (T2 + (B· _

= ui=1 '

2:7:1 (T2

'

J.l

)2) --1

+ (Bi- J.L)2)-1

as our Monte Carlo estimate of E,.(Bix). Note that (7.8) and (7.9) are applied separately to both the numerator and denominator, but using the same sample of B's. It is unwise to assume that the problem has been completely solved. The sample of B's generated from N(x, CT 2 ) will tend to concentrate around x,

7.3 Monte Carlo Sampling

213

whereas to satisfactorily account for the contribution of the Cauchy prior to the posterior mean, a significant portion of the ()'s should come from the tails of the posterior distribution. It may therefore appear that it is perhaps better to express the posterior mean in the form

E,.(()ix)

=

J~

() exp (- ( 8;~)

2

) n( ()) d()

a

oo

f~oo exp (- ( 8~~)2)

'

n( ()) d()

then sample ()'s from Cauchy(JL, T) and use the approximation

'"'m () L..i=l

'"'m

L..i=l

i

(

exp -

exp

(

(8i -x)

2

~

)

(8i-x)2)

-~

·

However, this is also not totally satisfactory because the tails of the posterior distribution are not as heavy as those of the Cauchy prior, and hence there will be excess sampling from the tails relative to the center. The implication is that the convergence of the approximation is slower and hence a larger error in approximation (for a fixed m). Ideally, therefore, sampling should be from the posterior distribution itself for a satisfactory approximation. With this view in mind, a variation of the above theme has been developed. This is called the Monte Carlo importance sampling. Consider (7.8) again. Suppose that it is difficult or expensive to sample directly from j, but there exists a probability density u that is very close to f from which it is easy to sample. Then we can rewrite (7.8) as

EJh(X)

=

L

h(x)f(x) dx

r h(x) u(x) f(x) u(x) dx

= Jx =

L{

=

Eu { h(X)w(X)} ,

h(x)w(x)} u(x) dx

where w(x) = f(x)ju(x). Now apply (7.9) with f replaced by u and h replaced by hw. In other words, generate i.i.d. observations X 1 , X 2 , ... from the density u and compute

The sampling density u is called the importance function. We illustrate importance sampling with the following example.

214

7 Bayesian Computations

Example 7.5. Suppose X 1,X2, ... ,Xn are i.i.d. N(B,a 2), where both Band a 2 are unknown. Independent priors are assumed for B and a 2 , where B has a double exponential distribution with density exp( -IBI)/2 and a 2 has the prior density of (1 + a 2)- 2. Neither of these is a standard prior, but robust choice of proper prior all the same. If the posterior mean of B is of interest, then it is necessary to compute

Because 1r(B, a 2lx) is not a standard density, let us look for a standard density close to it. Letting x denote the mean of the sample x 1, x 2 , ... , Xn and s;, = :2::::~= 1 (xi- x) 2 /n, note that

s;}) exp( -IBI)(1 + a 2)- 2 [s; + (B- x)2r/2+1 (a2)-Cn/2+2) exp (- 2;2 { (B- x)2 + s;}) x { [s; + (B- x)2r(n/2+1)} exp( -IBI)( a2 )2

1r(B, a 21x) ex (a 2)-nf 2 exp (- ; 2 { (B- x) 2 + 2 =

1:

2

a2

ex u1 (a 2IB)u2(B) exp( -IBI)(--2 ) 2 , 1+a where u 1(a 2 IB) is the density of inverse Gamma with shape parameter n/2 + 1 and scale parameter ~{(B- x) 2 + s;,}, and u 2 is the Student's t density with d.f. n + 1, location x and scale a multiple of sn. It may be noted that the tails of exp( -IBI)( 1 ~: 2 ) 2 do not have much of an influence in the presence of u 1(a 2 IB)u 2(B). Therefore, u( B, a 2) = u 1(a 2IB)u 2(B) may be chosen as a suitable importance function. This involves sampling B first from the density u 2 (B), and given this B, sampling a 2 from u 1(a 2 IB). This is repeated to generate further values of (B, a 2 ). Finally, after generating m of these pairs (B, a 2 ), the required posterior mean of B is approximated by

In some high-dimensional problems, a combination of numerical integration, Laplace approximation and Monte Carlo sampling seems to give appealing results. Delampady et al. (1993) use a Laplace-type approximation to obtain a suitable importance function in a high-dimensional problem. One area that we have not touched upon is how to generate random deviates from a given probability distribution. Clearly, this is a very important subject being the basis of any Monte Carlo sampling technique. Instead of providing a sketchy discussion from this vast area, we refer the reader to

7.4 Markov Chain Monte Carlo Methods

215

the excellent book by Robert and Casella (1999). We would, however, like to mention one recent and very important development in this area. This is the discovery of a very efficient algorithm to generate a sequence of uniform random deviates with a very big period of 2 19937 - 1. This algorithm, known as the Mersenne twister (MT), has many other desirable features as well, details on which may be found in Matsumoto and Nishimura (1998). The property of having a very large period is especially important because Monte Carlo simulation methods, especially MCMC, require very long sequences of random deviates for proper implementation.

7.4 Markov Chain Monte Carlo Methods 7.4.1 Introduction A severe drawback of the standard Monte Carlo sampling or Monte Carlo importance sampling is that complete determination of the functional form of the posterior density is needed for their implementation. Situations where posterior distributions are incompletely specified or are specified indirectly cannot be handled. One such instance is where the joint posterior distribution of the vector of parameters is specified in terms of several conditional and marginal distributions, but not directly. This actually covers a very large range of Bayesian analysis because a lot of Bayesian modeling is hierarchical so that the joint posterior is difficult to calculate but the conditional posteriors given parameters at different levels of hierarchy are easier to write down (and hence sample from). For instance, consider the normal-Cauchy problem of Example 7.1. As shown later in Section 7.4.6, this problem can be given a hierarchical structure wherein we have the normal model, the conjugate normal prior in the first stage with a hyperparameter for its variance and this hyperparameter again has the conjugate prior. Similarly, consider Example 7.2 where we have independent observations Xi ""' Poisson(Bi)· Now suppose the prior on the 8/s is a conjugate mixture. We again see (Problem 14) that a hierarchical prior structure can lead to analytically tractable conditional posteriors. It turns out that it is indeed possible in such cases to adopt an iterative Monte Carlo sampling scheme, which at the point of convergence will guarantee a random draw from the target joint posterior distribution. These iterative Monte Carlo procedures typically generate a random sequence with the Markov property such that this Markov chain is ergodic with the limiting distribution being the target posterior distribution. There is actually a whole class of such iterative procedures collectively called Markov chain Monte Carlo (MCMC) procedures. Different procedures from this class are suitable for different situations. As mentioned above, convergence of a random sequence with the Markov property is being utilized in this procedure, and hence some basic understanding of Markov chains is required. This material is presented below. This

216

7 Bayesian Computations

discussion as well as the following sections are mainly based on Athreya et al. (2003).

7.4.2 Markov Chains in MCMC A sequence of random variables {Xn}n;:=:o is a Markov chain if for any n, given the current value, Xn, the past {Xj, j :::; n -1} and the future { Xj : j ;::: n + 1} are independent. In other words,

P(A n BJXn) = P(AJXn)P(BJXn),

(7.10)

where A and B are events defined respectively in terms of the past and the future. Among Markov chains there is a subclass that has wide applicability. They are Markov chains with time hom*ogeneous or stationary transition probabilities, meaning that the probability distribution of Xn+ 1 given Xn = x, and the past, Xj : j :::; n- 1 depends only on x and does not depend on the values of Xj : j :::; n- 1 or n. If the set S of values {Xn} can take, known as the state space, is countable, this reduces to specifying the transition probability matrix P = ((pij)) where for any two values i, j in S, Pij is the probability that Xn+ 1 = j given Xn = i, i.e., of moving from state ito state j in one time unit. For state space S that is not countable, one has to specify a transition kernel or transition function P(x, ·) where P(x, A) is the probability of moving from x into A in one step, i.e., P(Xn+ 1 E AJXn = x). Given the transition probability and the probability distribution of the initial value X 0 , one can construct the joint probability distribution of { Xj : 0 :::; j :::; n} for any finite n. For example, in the countable state space case

P(Xo = io, X1 = i1, ... , Xn-1 = in-1, Xn =in) = P(Xn = inJXo = io, ... , Xn-1 = in-d xP(Xo = io,X1 = i1, .. . Xn-1 = in-1)

= Pin~linP(Xo = io, · · ·, Xn-1 = in-1) = P(Xo = io)PioitPi 1 iz · · ·Pin~lin·

A probability distribution 1T is called stationary or invariant for a transition probability P or the associated Markov chain {Xn} if it is the case that when the probability distribution of X 0 is 1T then the same is true for Xn for all n ;::: 1. Thus in the countable state space case a probability distribution 1T = {1Ti : i E S} is stationary for a transition probability matrix P if for each j inS, P(X 1 = j)

=L

P(X 1 = jJXo = i)P(Xo

= L 1TiPij = P(Xo

= j)

= i)

= 1Tj·

(7.11)

7.4 Markov Chain Monte Carlo Methods

In vector notation it says P with eigenvalue 1 and

7r =

( 1r 1 , 1r2 , ... )

is a left eigenvector of the matrix

7r=7rP.

(7.12)

Similarly, if S is a continuum, a probability distribution is stationary for the transition kernel P(-, ·) if

1r(A) =

217

L =is p(x) dx

1r

with density p(x)

P(x, A)p(x) dx

for all A C S. A Markov chain {Xn} with a countable state spaceS and transition probability matrix P (Pij)) is said to be irreducible if for any two states i and J. the probability of the Markov chain visiting j starting from i is positive, i.e., for some n 2: 1,p~7) P(Xn = jiXo = i) > 0. A similar notion of irreducibility, known as Harris or Doeblin irreducibility exists for the general state space case also. For details on this somewhat advanced notion as well as other results that we state here without proof, see Robert and Casella (1999) or Meyn and Tweedie (1993). In addition, Tierney(1994) and Athreya et al. (1996) may be used as more advanced references on irreducibility and MCMC. In particular, the last reference uses the fact that there is a stationary distribution of the Markov chain, namely, the joint posterior, and thus provides better and very explicit conditions for the MCMC to converge.

=( =

Theorem 7.6. (law of large numbers for Markov chains) Let {Xn}n::o:o be a Markov chain with a countable state space S and a transition probability matrix P. Further, suppose it is irreducible and has a stationary probability distribution 1r = (1ri : i E S) as defined in (7.11). Then, for any bounded function h : S--+ R and for any initial distribution of X 0 n-l

.!_ n

L

i=O

h(Xi) --+

L

h(j)Jrj

(7.13)

j

in probability as n--+ oo. A similar law of large numbers (LLN) holds when the state space S is not countable. The limit value in (7.13) will be the integral of h with respect to the stationary distribution 1r. A sufficient condition for the validity of this LLN is that the Markov chain { Xn} be Harris irreducible and have a stationary distribution Jr. To see how this is useful to us, consider the following. Given a probability distribution 1r on a set S, and a function h on S, suppose it is desired to compute the "integral of h with respect to 1r", which reduces to 2: 1 h(j)1r1 in the countable case. Look for an irreducible Markov chain {Xn} with state space S and stationary distribution 1r. Then, starting from some initial value

218

7 Bayesian Computations

X 0 , run the Markov chain {Xj} for a period of time, say 0, 1, 2, ... n - 1 and consider as an estimate

(7.14) By the LLN (7.13), this estimate fJn will be close to Lj h(j)nj for large n. This technique is called Markov chain Monte Carlo (MCMC). For example, if one is interested in n(A) = LjEA 'Trj for some A C S then by LLN (7.13) this reduces to 7rn(A)

1

n-1

=- L n o

IA(Xj)----+ n(A)

in probability as n----+ oo, where IA(Xj) = 1 if Xj E A and 0 otherwise. An irreducible Markov chain {Xn} with a countable state spaceS is called aperiodic if for some i E S the greatest common divisor, g.c.d. {n : p~~) > 0} = 1. Then, in addition to the LLN (7.13), the following result on the convergence of P(Xn = j) holds.

L IP(Xn = j)- 'Trjl----+ 0

(7.15)

j

as n ----+ oo, for any initial distribution of X 0 . In other words, for large n the probability distribution of Xn will be close to 1r. There exists a result similar to (7.15) for the general state space case that asserts that under suitable conditions, the probability distribution of Xn will be close to 1r as n----+ oo. This suggests that instead of doing one run of length n, one could do N independent runs each of length m so that n = N m and then from the ith run use only the mth observation, say, Xm,i and consider the estimate 1 N

MN,m

=N

L h(Xm,i)·

(7.16)

i=l

Other variations exist as well. Some of the special Markov chains used in MCMC are discussed in the next two sections.

7.4.3 Metropolis-Hastings Algorithm In this section, we discuss a very general MCMC method with wide applications. It will soon become clear why this important discovery has led to very considerable progress in simulation-based inference, particularly in Bayesian analysis. The idea here is not to directly simulate from the given target density (which may be computationally very difficult) at all, but to simulate an easy Markov chain that has this target density as the density of its stationary distribution. We begin with a somewhat abstract setting but very soon will get to practical implementation.

7.4 Markov Chain Monte Carlo Methods

219

Let S be a finite or countable set. Let 1r be a probability distribution on S. We shall call1r the target distribution. (There is room for slight confusion here because in our applications the target distribution will always be the posterior distribution, so let us note that 1r here does not denote the prior distribution, but just a standard notation for a generic target.) Let Q = ((Qij)) be a transition probability matrix such that for each i, it is computationally easy to generate a sample from the distribution {Qij : j E S}. Let us generate a Markov chain { Xn} as follows. If Xn = i, first sample from the distribution {qij : j E S} and denote that observation Yn. Then, choose Xn+l from the two values Xn and Yn according to

P(Xn+l P(Xn+l

=

=

YniXn, Yn)

XniXn, Yn)

=

=

p(Xn, Yn)

1- p(Xn, Yn),

(7.17)

where the "acceptance probability" p(·, ·) is given by . { 1fj Qji, pz,J ( . ") =mm

1}

(7.18)

1fi Qij

for all (i,j) such that 1fiQij > 0. Note that {Xn} is a Markov chain with transition probability matrix P = ( (Pij)) given by Qi]Pij Pij = {

1-

L Pik,

j # i, j = i.

(7.19)

kf'i

Q is called the "proposal transition probability" and p the "acceptance probability". A significant feature of this transition mechanism P is that P and satisfy HiPij = 1fjPji

for all i,j.

1r

(7.20)

This implies that for any j

(7.21) is a stationary probability distribution for P. Now assume that S is irreducible with respect to Q and Hi > 0 for all i in S. It can then be shown that Pis irreducible, and because it has a stationary distribution 1r, LLN (7.13) is available. This algorithm is thus a very flexible and useful one. The choice of Q is subject only to the condition that S is irreducible with respect to Q. Clearly, it is no loss of generality to assume that Hi > 0 for all i in S. A sufficient condition for the aperiodicity of P is that Pii > 0 for some i or equivalently

or,

1r

L #1

QijPij

< 1.

220

7 Bayesian Computations

A sufficient condition for this is that there exists a pair (i,j) such that 1ri% > 0 and 7rjqji < 1riqij. Recall that if P is aperiodic, then both the LLN (7.13) and (7.15) hold. If S is not finite or countable but is a continuum and the target distribution n( ·) has a density p( ·), then one proceeds as follows: Let Q be a transition function such that for each x, Q(x, ·) has a density q(x, y). Then proceed as in the discrete case but set the "acceptance probability" p(x, y) to be

. {p(y)q(y, x) } p(x, y) = mm p(x)q(x, y), 1 for all (x,y) such that p(x)q(x,y) > 0. A particularly useful feature of the above algorithm is that it is enough to know p( ·) up to a multiplicative constant as in the definition of the "acceptance probability" p(·, ·), only the ratios p(y) / p( x) need to be calculated. (In the discrete case, it is enough to know {ni} upto a multiplicative constant because the "acceptance probability" p(·, ·) needs only the ratios ndnj.) This assures us that in Bayesian applications it is not necessary to have the normalizing constant of the posterior density available for computation of the posterior quantities of interest.

7.4.4 Gibbs Sampling As was pointed out in Chapter 2, most of the new problems that Bayesians are asked to solve are high-dimensional. Applications to areas such as micro-arrays and image processing are some examples. Bayesian analysis of such problems invariably involve target (posterior) distributions that are high-dimensional multivariate distributions. In image processing, for example, typically one has N x N square grid of pixels with N = 256 and each pixel has k 2': 2 possible values. Thus each configuration has (256) 2 components and the state spaceS 2 has k( 256 ) configurations. To simulate a random configuration from a target distribution over such a large S is not an easy task. The Gibbs sampler is a technique especially suitable for generating an irreducible aperiodic Markov chain that has as its stationary distribution a target distribution in a highdimensional space but having some special structure. The most interesting aspect of this technique is that to run this Markov chain, it suffices to generate observations from univariate distributions. The Gibbs sampler in the context of a bivariate probability distribution can be described as follows. Let 1r be a target probability distribution of a bivariate random vector (X,Y). For each x, let P(x,·) be the conditional probability distribution of Y given X = x. Similarly, let Q(y, ·) be the conditional probability distribution of X given Y = y. Note that for each x, P(x, ·) is a univariate distribution, and for each y, Q(y, ·) is also a univariate distribution. Now generate a bivariate Markov chain Zn = (Xn, Yn) as follows: Start with some X 0 = x 0 . Generate an observation Yo from the distribution P(x 0 , ·). Then generate an observation X 1 from Q(Y0 , ·). Next generate an

7.4 Markov Chain Monte Carlo Methods

221

observation Y1 from P( X 1, ·) and so on. At stage n if Zn = (Xn, Yn) is known, then generate Xn+1 from Q(Yn, ·) and Yn+1 from P(Xn+1, ·). If 1r is a discrete distribution concentrated on {(xi, YJ) : 1 :S: i :S: K, 1 :S: j :S: L} and if 1fij = 1r(xi, YJ) then P(xi, YJ) = Hij /Hi. and

where 1fi. = Lj 1fij, 1f.j = Li 1fij· Thus the transition probability matrix R = ((r(ij),(k£))) for the {Zn} chain is given by '(ij),(k£) = Q(yj, xk)P(xk, Y£) 1fkj 1fk£ 1f.j 1fk.

It can be verified that this chain is irreducible, aperiodic, and has 1r as its stationary distribution. Thus LLN (7.13) and (7.15) hold in this case. Thus for large n, Zn can be viewed as a sample from a distribution that is close to 7f and one can approximate Li,j h(i,j)Hij by 2.:::7= 1 h(Xi, Y;)jn. As an illustration, consider sampling from (

~) ~

N2 ( (

i]).

~), [ ~

Note that the conditional distribution of X given Y = y and that of Y given X= x are X[Y = y ~ N(py, 1- p 2 ) and Y[X = x ~ N(px, 1- p 2 ).

(7.22)

Using this property, Gibbs sampling proceeds as described below to generate (Xn, Yn), n = 0, 1, 2, ... , by starting from an arbitrary value xo for Xo, and repeating the following steps fori = 0, 1, ... , n. 1. Given Xi for X, draw a random deviate from N(pxi, 1- p 2 ) and denote it by Y;. 2. Given Yi for Y, draw a random deviate from N(pyi, 1- p 2 ) and denote it by xi+1· The theory of Gibbs sampling tells us that if n is large, then (xn, Yn) is a random draw from a distribution that is close to N 2 (

(

~) , [ ~

i] ).

To see

why Gibbs sampler works here, recall that a sufficient condition for the LLN (7.13) and the limit result (7.15) is that an appropriate irreducibility condition holds and a stationary distribution exists. From steps 1 and 2 above and using (7.22), one has and xi+1

= pY; + .Jl=P2 ~i,

where 'T/i and ~i are independent standard normal random variables independent of Xi. Thus the sequence {Xi} satisfies the stochastic difference equation

222

7 Bayesian Computations

where ui+1 =

pyf1 - p 2 TJi

+ yi1- p2 ~i·

Because 7Ji,~i are independent N(O, 1) random variables, Ui+ 1 is also a normally distributed random variable with mean 0 and variance p 2 (1- p 2 ) + (1p2) = 1- p4 . Also {Ui}i2 1 being i.i.d., makes {Xi}i2 0 a Markov chain. It turns out that the irreducibility condition holds here. Turning to stationarity, note that if X 0 is a N(O, 1) random variable, then X 1 = p 2 X 0 + U1 is also a N(O, 1) random variable, because the variance of X 1 = p4 + 1- p4 = 1 and the mean of X 1 is 0. This makes the standard N(O, 1) distribution a stationary distribution for {Xn}· The multivariate extension of the above-mentioned bivariate case is very straightforward. Suppose 1r is a probability distribution of a k-dimensional random vector (X 1,Xz, ... ,Xk)· Ifu = (u 1,uz, ... ,uk) is any k-vector, let u_i = ( u 1, u 2, ... , Ui- 1, ui+ 1, ... , uk) be the k -1 dimensional vector resulting by dropping the ith component ui. Let ni(·lx-i) denote the univariate conditional distribution of Xi given that X_i (X 1,Xz,Xi_ 1,Xi+ 1, ... ,Xk) = X-i· Now starting with some initial value for X 0 = (x 01 , x 02 , ... , Xok) generate X 1 = (Xll,X 12 , ... ,X1k) sequentially by generating Xll according to the univariate distribution 7r1(·lxo_,) and then generating x12 according to 7rz(·I(Xll,xo3,Xo4,·· .,Xok) and so on. The most important feature to recognize here is that all the univariate conditional distributions, XiiX-i = x_i, known as full conditionals should easily allow sampling from them. This turns out to be the case in most hierarchical Bayes problems. Thus, the Gibbs sampler is particularly well adapted for Bayesian computations with hierarchical priors. This was the motivation for some vigorous initial development of Gibbs sampling as can be seen in Gelfand and Smith (1990). The Gibbs sampler can be justified without showing that it is a special case of the Metropolis-Hastings algorithm. Even if it is considered a special case, it still has special features that need recognition. One such feature is that full conditionals have sufficient information to uniquely determine a multivariate joint distribution. This is the famous Hammersley- Clifford theorem. The following condition introduced by Besag (1974) is needed to state this result.

=

Definition 7. 7. Let p(y 1, ... , Yk) be the joint density of a random vector Y = (Y1 , ... , Yk) and let p(i) (yi) denote the marginal density of Y;, i = 1, ... , k. If p(i) (Yi) > 0 for every i = 1, ... , k implies that p(y 1, . .. , Yk) > 0, then the joint density p is said to satisfy the positivity condition. Let us use the notation Pi(YiiY1, ... , Yi-1, Yi+1, ... , Yk) for the conditional density of Y;IY - i = y -i·

Theorem 7.8. (Hammersley-Clifford) Under the positivity condition, the joint density p satisfies

7.4 Markov Chain Monte Carlo Methods

223

for every y andy' in the support of p. Proof. For y and y' in the support of p, p(y1, · · · ,yk)

= Pk(YkiY1, · · · ,Yk-dP(Y1, · · · ,Yk-d _ Pk(YkiY1,···,Yk-d ( ') ('I )PY1,···,Yk-1,Yk Pk Yk Y1, · · ·, Yk-1

-

Pk(Yk IY1) ... ) Yk-d Pk-1 (Yk-11Y1, ... ) Yk-2, YU

Pk(Y~ IY1) ... ) Yk-d Pk-1 (Y~-11Y1, ... ) Yk-2, YU xp(y1, · · ·, Y~-1, YU

D

It can be shown also that under the positivity condition, the Gibbs sampler generates an irreducible Markov chain, thus providing the necessary convergence properties without recourse to the M-H algorithm. Additional conditions are, however, required to extend the above theorem to the non-positive case, details of which may be found in Robert and Casella (1999).

7.4.5 Rao-Blackwellization The variance reduction idea of the famous Rao-Blackwell theorem in the presence of auxiliary information can be used to provide improved estimators when MCMC procedures are adopted. Let us first recall this theorem. Theorem 7.9. (Rao-Blackwell theorem) Let o(X1 , X 2 , ... , Xn) be an estimator of with finite variance. Suppose that T is sufficient for and let o*(T), defined by o*(t) = E(o(X1 , X 2 , ... , Xn)IT = t), be the conditional expectation of o(X1 , X 2 , ... , Xn) given T = t. Then

e

The inequality is strict unless o ofT.

e,

= o*,

or equivalently, o is already a function

Proof. By the property of iterated conditional expectation,

Therefore, to compare the mean squared errors (MSE) of the two estimators, we need to compare their variances only. Now,

224

7 Bayesian Computations

Var(8(X 1 , X2, ... , Xn)) = Var [E(8IT)] = Var(8*)

+ E [Var(8IT)]

+ E [Var(8IT)]

> Var(8*), unless Var(8IT) = 0, which is the case only if 8 itself is a function ofT.

D

The Rao-Blackwell theorem involves two key steps: variance reduction by conditioning and conditioning by a sufficient statistic. The first step is based on the analysis of variance formula: For any two random variables SandT, because Var(S) = Var(E(SIT)) + E(Var(SIT)), one can reduce the variance of a random variable S by taking conditional expectation given some auxiliary information T. This can be exploited in MCMC. Let (Xj,Yj),j = 1,2, ... ,N be the data generated by a single run of the Gibbs sampler algorithm with a target distribution of a bivariate random vector (X, Y). Let h(X) be a function of the X component of (X, Y) and let its mean value be J-l· Suppose the goal is to estimate J-l. A first estimate is the sample mean of the h(Xj),j = 1, 2, ... , N. From the MCMC theory, it can be shown that as N -+ oo, this estimate will converge to J-l in probability. The computation of variance of this estimator is not easy due to the (Markovian) dependence of the sequence {Xj,j = 1, 2, ... , N}. Now suppose we make n independent runs of Gibbs sampler and generate (Xij, Yij ), j = 1, 2, ... , N; i = 1, 2, ... , n. Now suppose that N is sufficiently large so that (X iN, YiN) can be regarded as a sample from the limiting target distribution of the Gibbs sampling scheme. Thus (XiN, YiN), i = 1, 2, ... , n are i.i.d. and hence form a random sample from the target distribution. Then one can offer a second estimate of J-t-the sample mean of h(XiN ), i = 1, 2, ... , n. This estimator ignores a good part of the MCMC data but has the advantage that the variables h(XiN ), i = 1, 2, ... , n are independent and hence the variance of their mean is of order n- 1 . Now applying the variance reduction idea of the Rao-Blackwell theorem by using the auxiliary information YiN, i = 1, 2, ... , n, one can improve this estimator as follows: Let k(y) = E(h(X)IY = y). Then for each i, k(YiN) has a smaller variance than h(XiN) and hence the following third estimator, 1 n - Lk(YiN), n i=l

has a smaller variance than the second one. A crucial fact to keep in mind here is that the exact functional form of k(y) be available for implementing this improvement.

7.4 Markov Chain Monte Carlo Methods

225

7.4.6 Examples

Example 7.10. (Example 7.1 continued.) Recall that XIB"' N(B, 0' 2 ) with known 0' 2 and e "' Cauchy(f.l, T). The task is to simulate e from the posterior distribution, but we have already noted that sampling directly from the posterior distribution is difficult. What facilitates Gibbs sampling here is the result that the Student's t density, of which Cauchy is a special case, is a scale mixture of normal densities, with the scale parameter having a Gamma distribution (see Section 2. 7.2, Jeffreys test). Specifically,

so that 7r( e) may be considered the marginal prior density from the joint prior density of (B, A) where BIA"' N(ft, T 2 /A) and A"' Gamma(1/2, 1/2).

It can be noted that this leads to an implicit hierarchical prior structure with A being the hyperparameter. Consequently, 1r(Bix) may be treated as the marginal density from 1r(B, Alx). Now note that the full conditionals of 1r(B, Alx) are standard distributions from which sampling is easy. In particular, (7.23) (7.24) Thus, the Gibbs sampler will use (7.23) and (7.24) to generate (B, A) from

1r(B, Alx). Example 7.11. Consider the following example due to Casella and George given in Arnold (1993). Suppose we are studying the distribution of the number of defectives X in the daily production of a product. Consider the model (X I Y, e)"' binomial(Y, B), where Y, a day's production, is a random variable with a Poisson distribution with known mean A, and is the probability that any product is defective. The difficulty, however, is that Y is not observable, and inference has to be made on the basis of X only. The prior distribution is such that (B I Y = y) "' Beta( a, 1), with known a and 1 independent of Y. Bayesian analysis here is not a particularly difficult problem because the posterior distribution of BIX = x can be obtained as follows. First, note that XIB "'Poisson(AB). Next, Beta( a, I)· Therefore,

e

e"'

(7.25)

226

7 Bayesian Computations

The only difficulty is that this is not a standard distribution, and hence posterior quantities cannot be obtained in closed form. Numerical integration is quite simple to perform with this density. However, Gibbs sampling provides an excellent alternative. Instead of focusing on OIX directly, view it as a marginal component of (Y, (} I X). It can be immediately checked that the full conditionals of this are given by YIX = x, (} "'x + Poisson(.X(1 - 0) ), and OIX = x, Y = y"' Beta( a+ x, 1 + y- x) both of which are standard distributions. Example 7.12. (Example 7.11 continued.) It is actually possible here to sample from the posterior distribution using what is known as the accept-reject Monte Carlo method. This widely applicable method operates as follows. Let g(x)/ K be the target density, where K is the possibly unknown normalizing constant of the unnormalized density g. Suppose h(x) is a density that can be simulated by a known method and is close to g, and suppose there exists a known constant c > 0 such that g(x) < ch(x) for all x. Then, to simulate from the target density, the following two steps suffice. (See Robert and Casella (1999) for details.) Step 1. Generate Y"' hand U"' U(O, 1); Step 2. Accept X= Y if U::; g(Y)/{ch(Y)}; return to Step 1 otherwise. The optimal choice for cis sup{g(x)/h(x), but even this choice may result in undesirably large number of rejections. In our example, from (7.25),

g(O) = exp(-.XO)Ox+o-l(1-

0)'~'- 1 1{0::;

(}::; 1},

so that h(O) may be chosen to be the density of Beta(x +a,,). Then, with the above-mentioned choice for c, if(}"' Beta(x+a,/) is generated in Step 1, its 'acceptance probability' in Step 2 is simply exp(-.XO). Even though this method can be employed here, we, however, would like to use this technique to illustrate the Metropolis-Hastings algorithm. The required Markov chain is generated by taking the transition density q(z, y) = q(ylz) = h(y), independently of z. Then the acceptance probability is

p(z, y)

=

. {g(y)h(z) } mm g(z)h(y), 1

=min {exp (-.X(y- z)), 1}. Thus the steps involved in this "independent" M-H algorithm are as follows. Start at t = 0 with a value x 0 in the support of the target distribution; in this case, 0 < x 0 < 1. Given Xt, generate the next value in the chain as given below. (a) Draw Yt from Beta(x +a,,). (b) Let x = { Yt with probability Pt (t+l) Xt otherwise,

7.4 Markov Chain Monte Carlo Methods

227

where Pt = min{ exp ( ->.(yt- Xt)), 1}. (c) Set t = t + 1 and go to step (a). Run this chain until t = n, a suitably chosen large integer. Details on its convergence as well as why independent M-H is more efficient than acceptreject Monte Carlo can be found in Robert and Casella (1999). In our example, for x = 1, a= 1, 'Y = 49 and>.= 100, we simulated such a Markov chain. The resulting frequency histogram is shown in Figure 7.1, with the true posterior density super-imposed on it. Example 7.13. In this example, we discuss the hierarchical Bayesian analysis of the usual one-way ANOVA. Consider the model Yij

=

()i

Eij

rv

N(O,O"l),j

+Eij,j = 1, ... ,ni;i = 1, ... ,k; = 1, ...

,ni;i

= 1, ...

,k,

(7.26)

and

0"7

be such that they

and are independent. Let the first stage prior on are i.i.d. with (}i

N(f.ktr, 0";),

rv

i

=

()i

1, ... , k;

(]"7 "'inverse Gamma(a 1 , bl),

i = 1, ... , k.

The second stage prior on f.ktr and O"; is

-

/ ::=r\1\

f

~

~

1\ -

1\c\

~

~~

"ltffl, 0

0.02

Ll

0.04

Fig. 7.1. M-H frequency histogram and true posterior density.

0.06

228

7 Bayesian Computations

Here a 1 , a 2 , b1 , b2 , J.Lo, and a5 are all specified constants. Let us concentrate on computing Sufficiency reduces this to considering only -

Yi

1

n;

""" Yij, i = 1, ... , k

= -

n~ t

and

j=l

n;

""" si2= ~(Yij-

-2· Yi) ,z =

1, ... , k.

j=l

From normal theory, 2 2 Yil8, u "'N(Oi, ai /ni),

i

=

1, ... , k,

which are independent and are also independent of STI8,u 2 i=1, ... ,k,

rvalx;;-ll

which again are independent. To utilize the Gibbs sampler, we need the full conditionals of7r(8,u 2 ,J.L7r,a;ly). It can be noted that it is sufficient, and in fact advantageous to consider the conditionals, (i) 7r(8lu 2 ,J.L7r,a;,y), (ii) 1r(u2 l8, fL1r, y), (iii) 7r(J.L7rla;,8,u 2 ,y), and (iv) 7r(a;IJ.L7r, 8, u 2 , y), rather than considering the set of all univariate full conditionals because of the special structure in this problem. First note that

a;,

and hence

(7.27) which determines (i). Next we note that, given 8, from (7.26), St 2 = I:,j~ 1 (Yij - Oi) 2 is sufficient for af, and they are independently distributed. Thus we have, and are independent, and

a; are i.i.d. inverse Gamma (a

1,

bl), so that

7.4 Markov Chain Monte Carlo Methods

229 (7.28)

and they are independent fori= 1, ... , k, which specifies (ii). Turning to the full conditional of Jl1r, we note from the hierarchical structure that the conditional distribution of Jl1r I u 2 , y is the same as the conditional distribution of Jl1r 10';, (). To determine this distribution, note that

0';, (),

ei IJL1r, 0'; "" N(JL1f, 0'; ), fori= 1, ... , k and are i.i.d. and Jl1r ""N(Jlo, 0'6). Therefore, treating() to be a random sample from N(J11f, so that 1J = ~7= 1 B;jk is sufficient for J11r, we have the joint distribution,

0';),

Thus we obtain,

which provides (iii). Just as in the previous case, the conditional distribution of IJ11r, (), u 2 , y turns out to be the same as the conditional distribution of IJ11r, (). To obtain this, note again that

0';

0';

for i = 1, 'k and are i.i.d. so that this time ~7=1 Further

0';.

k

L(Bi- Jl1r ) 210'; ""O';X~ and

(ei - jl1f ) 2 is sufficient for

0'; ""inverse Gamma(a2, b2),

i=l

so that

This gives us (iv), thus completing the specification of all the required full conditionals. It may be noted that the Gibbs sampler in this problem requires simulations from only the standard normal and the inverse Gamma distributions.

Reversible Jump MCMC There are situations, especially in model selection problems, where the MCMC procedure should be capable of moving between parameter spaces of different dimensions. The standard M-H algorithm described earlier is incapable

230

7 Bayesian Computations

of such movements, whereas the reversible jump algorithm of Green (1995) is an extension of the standard M-H algorithm to allow exactly this possibility. The basic idea behind this technique as applied to model selection is as follows. Given two models M 1 and M 2 with parameter sets fh and fh, which are possibly of different dimensions, fill the difference in the dimensions by supplementing the parameter sets of these models. In other words, find auxiliary variables 112 and 12 1 such that (8 1,1 12) and (82,12 1) can be mapped with a bijection. Now use the standard M-H algorithm to move between the two models; for moves of the M-H chain within a model, the auxiliary variables are not needed. We sketch this procedure below, but for further details refer to Robert and Casella (1999), Green (1995), Sorensen and Gianola (2002), Waagepetersen and Sorensen (2001), and Brooks et al. (2003). Consider models M 1 , M 2 , ... where model Mi has a continuous parameter space fh The parameter space for the model selection problem as a whole may be taken to be

Let f(xiMi, Oi) be the model density under model Mi, and the prior density be

where 1ri is the prior probability of model Mi and n(OiiMi) is the prior density conditional on Mi being true. Then the posterior probability of any B C U/7] is

where is the posterior density restricted to Mi. To compute the Bayes factor of Mk relative to Mz, we will need

P"(Mklx) nz P"(Mzlx) 7rk' where

P" ( Mi lx)

ni .fte. n(OiiMi)f(xiMi, Oi) dOi

= =-----'+'-~__,----,----,---------:---,------c-

l::j 7rj fej n( Oj IMj) f(xiMj, Oj) dOj

is the posterior probability of Mi. Therefore, for the target density n(Oix), we need a version of the M-H algorithm that will facilitate the above-shown computations. Suppose Oi is a vector of length ni. It suffices to focus on moves between Oi in model Mi and Oj in model Mj with ni < nj· The scheme provided by Green (1995) is as follows. If the current state of the chain is (Mi,Oi), a new value (Mj,Oj) is proposed for the chain from a proposal

7.4 Markov Chain Monte Carlo Methods

231

(transition) distribution Q(8i, d8j ), which is then accepted with a certain acceptance probability. To move from model Mi to Mj, generate a random vector V of length nj - ni from a proposal density nj-ni

IT

'1/Jij(v) =

'1/J(vrn).

rn=l

Identify an appropriate bijection map

and propose the move from 8i to 8j using 8j = hij(8i, V). The acceptance probability is then

where

. (8 8 ) a,J

"

1

-

7r(8jiMj,X)Pji(8j) Iahij(8i, v) I 7r(8iiMi,x)Pij(8i)'I/Jij(v) 8(8i, v) '

with Pij(8i) denoting the (user-specified) probability that a proposed jump to model Mj is attempted at any step starting from 8i E th Note that "L, j Pii = 1.

Example 1.14. For illustration purposes, consider the simple problem of comparing two normal means as in Sorensen and Gianola (2002). Then, the two models to be compared are

To implement the reversible jump M-H algorithm we need the map, h 12 taking (v, o- 2 , V) to (v1 , v 2 , o- 2 ). A reasonable choice for this is the linear map

7 .4. 7 Convergence Issues As we have already seen, Monte Carlo sampling based approaches for inference make use of limit theorems such as the law of large numbers and the central limit theorem to justify their validity. When we add a further dimension to this sampling and adopt MCMC schemes, stronger limit theorems are needed. Ergodic theorems for Markov chains such as those given in equations (7.13) and (7.15) are these useful results. It may appear at first that this procedure

232

7 Bayesian Computations

necessarily depends on waiting until the Markov chain converges to the target invariant distribution, and sampling from this distribution. In other words, one needs to start a large number of chains beginning with different starting points, and pick the draws after letting these chains run sufficiently long. This is certainly an option, but the law of large numbers for dependent chains, (7.13) says also that this is unnecessary, and one could just use a single long chain. It may, however, be a good idea to use many different chains to ensure that convergence indeed takes place. For details, see Robert and Casella (1999). There is one important situation, however, where MCMC sampling can lead to absurd inferences. This is where one resorts to MCMC sampling without realizing that the target posterior distribution is not a probability distribution, but an improper one. The following example is similar to the normal problem (see Exercise 13) with lack of identifiability of parameters shown in Carlin and Louis (1996).

Example 7.15. (Example 7.11 continued.) Recall that, in this problem, (X

I

Y, e) "' binomial(Y, e), where Y I A "' Poisson( A). Earlier, we worked with

a known mean A, but let us now see if it is possible to handle this problem with unknown A. Because Y is unobservable and only X is observable, there exists an 'identifiability' problem here, as can be seen by noting that Xle "' Poisson(,\e). We already have the Beta(a,/) prior on e. Suppose 0

e < 1,,\ > 0.

(7.31)

This joint density is improper because

11

1

00

=

=

=

exp(-Ae)Axex+a- 1(1-e)'- 1d,\de

11 (100 exp( -Ae)Ax d,\) ex+a-1(1- ep-1 de

t

lo

r(x + 1) ex+a-1(1- e)'-1 de ex+1

r(x

+ 1) 1

1

e"'- 2 (1- ep- 1 de

= 00.

In fact, the marginal distributions are also improper. However, it has full conditional distributions that are proper:

Ale, X"' Gamma (x

+ 1, e)

and 7r(eiA, x) ex exp( -Ae)ex+<>- 1(1- e)'- 1.

Thus, for example, the Gibbs sampler can be successfully employed with these proper full conditionals. To generate from 1r(el,\,x), one may use the independent M-H algorithm described in Example 7.12. Any inference on the

e

7.5 Exercises

233

marginal posterior distributions derived from this sample, however, will be totally erroneous, whereas inferences can indeed be made on >..e. In fact, the non-convergence of the chain encountered in the above example is far from being uncommon. Often when we have a hierarchical prior, the prior at the final stage of the hierarchy is an improper objective prior. Then it is not easy to check that the joint posterior is proper. Then none of the theorems on convergence of the chains may apply, but the chain may yet seem to converge. In such cases, inference based on MCMC may be misleading in the sense of what was seen in the example above.

7.5 Exercises 1. (Flury and Zoppe (2000)) A total of m + n lightbulbs are tested in two independent experiments. In the first experiment involving n lightbulbs, the exact lifetimes Y1, ... , Yn of all the bulbs are recorded. In the second involving m lightbulbs, the only information available is whether these lightbulbs were still burning at some fixed time t > 0. This is known as right-censoring. Assume that the distribution of lifetime is exponential with mean 1je, and use n(e) ex: e- 1 . Find the posterior mode using the E-M algorithm. 2. (Flury and Zoppe (2000)) In Problem 1, use uniform(O, e) instead of exponential for the lifetime distribution, and n(e) = I(o,oo)(e). Show that the E-M algorithm fails here if used to find the posterior mode. 3. (Inverse c.d.f. method) Show that, if the c.d.f. F(x) of a random variable X is continuous and strictly increasing, then U = F(X) "' U[O, 1], and if V "' U[O, 1], then Y F- 1 (V) has c.d.f F. Using this show that if U "' U[O, 1], -ln U / (3 is an exponential random variable with mean (3- 1 . 4. (Box-Muller transformation method) Let U1 and U2 be a pair of independent Uniform (0, 1) random variables. Consider first a transformation to

=

and then let X

= R cos V; Y = R sin V.

Show that X and Y are independent standard normal random variables. 5. Prove that the accept-reject Monte Carlo method given in Example 7.12 indeed generates samples from the target density. Further show that the expected number of draws required from the 'proposal density' per observation is c- 1 . 6. Using the methods indicated in Exercises 1, 2, and 3 above, or combinations thereof, prove that the standard continuous probability distributions can be simulated.

234

7 Bayesian Computations

7. Consider a discrete probability distribution that puts mass p; on point x;, i = 0, 1, .... Let U"" U(O, 1), and define a new random variable Y as follows. y _ { xo if U ::::; Po; . x; if 2:~:~Pj < U::::; L~=oPj, i ~ 1. What is the probability distribution of Y? 8. Show that the random sequence generated by the independent M-H algorithm is a Markov chain. 9. (Robert and Casella (1999)) Show that the Gamma distribution with a non-integer shape parameter can be simulated using the accept-reject method or the independent M-H algorithm. 10. Gibbs Sampling for Multinomial. Consider the ABO Blood Group problem from Rao (1973). The observed counts in the four blood groups, 0, A, B, and ABare as given in Table 7.3. Assuming that the inheritance of these blood groups is controlled by three alleles, A, B, and 0, of which 0 is recessive to A and B, there are six genotypes 00, AO, AA, BO, BB, and AB, but only four phenotypes. If r, p, and q are the gene frequencies of 0, A, and B, respectively (with p + q + r = 1), then the probabilities of the four phenotypes assuming Hardy-Weinberg equilibrium are also as shown in Table 7.3. Thus we have here a 4-cell multinomial probability vector that is a function of three parameters p, q, r with p + q + r = 1. One may wish to formulate a Dirichlet prior for p, q, r. But it will not be conjugate to the 4-cell multinomial likelihood function in terms of p, q, r from the data, and this makes it difficult to work out the posterior distribution of p, q, r. Although no data are missing in the real sense of the term, it is profitable to split each of the nA and nB cells into two: nA into nAA, nAo with corresponding probabilities p 2 , 2pr and nB into nEB, nBo with corresponding probabilities q2 , 2qr, and consider the 6-cell multinomial problem as a complete problem with nAA, nEB as 'missing' data. Table 7.3. ABO Blood Group Data

Cell Count Probability no= 176 nA 182

ns

= =

nAB=

60 17

r

+ 2pr l + 2qr

p

2

2pq

= no+ nA + nB +nAB, and denote the observed data by n = (no, nA, nB, nAB). Consider estimation of p, q, r using a Dirichlet prior with parameters a, {3, '"Y with the 'incomplete' observed data n. Let N

The likelihood upto a multiplicative constant is

L(p, q, r)

=

r2no (p2 + 2pr)nA (q2 + 2qr)nB (pq)nAB.

7.5 Exercises

235

The posterior density of (p, q, r) given n is proportional to

r2no+!-l(p2

+ 2pr)nA(q2 + 2qr)nB(p)nAB+a-l(q)nAB+.6-1.

Let nA = nAA + nAo, ns = nss + nso, and write noo for no. Verify that if we have the 'complete' data, ii = (noo, nAA, nAo, nss, nso, nAs), then the likelihood is, upto a multiplicative constant

where

n1 = nAA

1

1

+ 2nAB + 2nAo

+ nB

= 2nAB + nss + 2nso

n6

= 2nAo + 2nso + noo.

1

1

1

1

Show that the posterior distribution of (p, q, r) given ii is Dirichlet with parameters n1 +a- 1, nt + f3- 1, n6 + 1- 1, when the prior is Dirichlet with parameters (a, /3, 1). Show that the conditional distributions of (nAA, nss) given ii and (p, q, r) is that of two independent binomials: 2

(nAAin,p, q, r)

rv

binomial(nA,

p2

p ), +2pr

2

(nssln,p,q,r) (p, q, rln, nAA, nss)

rv

rv

binomial(ns,

q ), and q +2qr 2

Dirichlet(n1 +a- 1, nt

+ /3-

1, n6

+ 1- 1).

Show that the Rao-Blackwellized estimate of (p, q, r) from a Gibbs sample of size m is

where the superscript i denotes the ith draw. 11. (M-H for the Weibull Model: Robert (2001)). The following twelve observations are from a simulated reliability study: 0.56, 2.26, 1.90, 0.94, 1.40, 1.39, 1.00, 1.45, 2.32, 2.08, 0.89, 1.68. A Wei bull model with the following density form is considered appropriate:

f(xla, ry) ex aryx"'- 1e-ryx"', 0 < x < oo, with parameters (a,ry). Consider the prior distribution

236

7 Bayesian Computations

7r( a, 17) ex e -a17{3-l e -t;1). The posterior distribution of (a,ry) given the data (x 1 ,x 2 , ... ,xn) has density

To get a sample from the posterior density, one may use the M-H algorithm with proposal density

1 {-a'a- - -11'17 } ,

q(a', ry'la, 17) = - exp

ary

which is a product of two independent exponential distributions with means a, 11· Compute the acceptance probability p((a',ry'), (o:(tl,ry(tl)) at the tth step of the M-H chain, and explain how the chain is to be generated. 12. Complete the construction of the reversible jump M-H algorithm in Example 7.14. In particular, choose an appropriate prior distribution, proposal distribution and compute the acceptance probabilities. 13. (Carlin and Louis (1996)) Suppose Yl,Y2, ... ,Yn is an i.i.d. sample with

Yil81, 82 '""N(81 + 82, a 2), where a 2 is assumed to be known. Independent improper uniform prior distributions are assumed for 81 and 82. (a) Show that the posterior density of (81,82IY) is

7r(8l' 82IY) ex exp( -n(8l + 82 - y) 2I (2a 2) )I( (81' 82) E R 2), which is improper, integrating to oo (over R 2 )). (b) Show that the marginal posterior distributions are also improper. (c) Show that the full conditional distributions of this posterior distribution are proper. (d) Explain why a sample generated using the Gibbs sampler based on these proper full conditionals will be totally useless for any inference on the marginal posterior distributions, whereas inferences can indeed be made on 81 +82. 14. Suppose X 1 , X 2 , ... , Xk are independent Poisson counts with Xi having mean 8i. 8i are a priori considered related, but exchangeable, and the prior k

7r(81, ... '8k) ex (1 +

L 8i)-(k+l)'

i=l is assumed. (a) Show that the prior is a conjugate mixture. (b) Show how the Gibbs sampler can be employed for inference.

7.5 Exercises

237

15. Suppose X 1 ,X2, ... ,Xn are i.i.d. random variables with

and independent scale-invariant non-informative priors on >. 1 and >. 2 are used. i.e., 1r(>.1, >-2) ex (>.1>-2)- 1I(>-1 > 0, A2 > 0). (a) Show that the marginals of the posterior, 1r(>. 1, >- 21x) are improper, but the full conditionals are standard distributions. (b) What posterior inferences are possible based on a sample generated from the Gibbs sampler using these full conditionals?

8 Some Common Problems in Inference

We have already discussed some basic inference problems in the previous chapters. These include the problems involving the normal mean and the binomial proportion. Some other usually encountered problems are discussed in what follows.

8.1 Comparing Two Normal Means Investigating the difference between two mean values or two proportions is a frequently encountered problem. Examples include agricultural experiments where two different varieties of seeds or fertilizers are employed, or clinical trials involving two different treatments. Comparison of two binomial proportions was considered in Example 4.6 and Problem 8 in Chapter 4. Comparison of two normal means is discussed below. Suppose the model for the available data is as follows. Yn, ... , Y1 n, is a random sample of size n 1 from a normal population, N(B 1,err), whereas Y21, ... , Y2n 2 is an independent random sample of size n2 from another normal population, N(B 2 , er§). All the four parameters 8 1, 8 2, err, and er§ are unknown, but the quantity of inferential interest is 7) = el - e2. It is convenient to consider the case, err = er§ = er 2 separately. In this case, (Yi, Y2,s 2) is jointly sufficient for (B1,B2,er 2) where s 2 = C2=~~ 1 (Yli- Y1? + 2::::7~ 1 (Y2j- Y2) 2)/(n1 + n2- 2). Further, given (8 1, 82, er 2),

and they are independently distributed. Upon utilizing the objective prior, n(B 1,B2,er 2) ex er- 2, one obtains

and hence

240

8 Some Common Problems in Inference

(8.1) Now, note that

Consequently, integrating out a 2 from (8.1) yields,

(8.2) or, equivalently

ar

In many situations, the assumption that = a~ is not tenable. For example, in a clinical trial the populations corresponding with two different treatments may have very different spread. This problem of comparing means when we have unequal and unknown variances is known as the Behrens-Fisher problem, and a frequentist approach to this problem has already been discussed in Problem 17, Chapter 2. We discuss the Bayesian approach now. We have that (l\' si) is sufficient for (el' a?) and (Y2' s~) is sufficient for (e2' a~)' where s7 = I:7~ 1 (Yi1- fi) 2/(ni- 1), i = 1, 2. Also, given (0 1 , 02, af,

an,

and further, they are all independently distributed. Now employ the objective prior and proceed exactly as in the previous case. It then follows that under the posterior distribution also el and e2 are independent, and that and It may be immediately noted that the posterior distribution of TJ = el - e2, however, is not a standard distribution. Posterior computations are still quite easy to perform because Monte Carlo sampling is totally straightforward. Simply generate independent deviates el and e2 repeatedly from (8.3) and utilize the corresponding TJ = el - e2 values to investigate its posterior distribution. Problem 4 is expected to apply these results. Extension to the k-mean problem or one-way ANOVA is straightforward. A hierarchical Bayes approach to this problem and implementation using MCMC have already been discussed in Example 7.13 in Chapter 7.

8.2 Linear Regression

241

8.2 Linear Regression We encountered normal linear regression in Section 5.4 where we discussed prior elicitation issues in the context of the problem of inference on a response variable Y conditional on some predictor variable X. Normal linear models in general, and regression models in particular are very widely used. We have already seen an illustration of this in Example 7.13 in Section 7.4.6. We intend to cover some of the important inference problems related to regression in this section. Extending the simple linear regression model where E(YI,Bo, ,61 , X = x) = ,60 + ,61 x to the multiple linear regression case, E(Yif3, X = x) = ,6'x, yields the linear model (8.4) y=Xf3+E, where y is the n-vector of observations, X the n x p matrix having the appropriate readings from the predictors, ,13 the p-vector of unknown regression coefficients, and E the n-vector of random errors with mean 0 and constant variance a- 2 . The parameter vector then is (,13, a- 2 ), and most often the statistical inference problem involves estimation of ,13 and also testing hypotheses involving the same parameter vector. For convenience, we assume that X has full column rank p < n. We also assume that the first column of X is the vector of 1 's, so that the first element of ,13, namely ,61 , is the intercept. If we assume that the random errors are independent normals, we obtain the likelihood function for (,13, a- 2 ) as

f(ylf3, a-

2

)

=

[~a- r

exp {-

=[~o-r exp{- 2 ~ 2

2 ~ 2 (y- Xj3)'(y- X,6)}

[(y-y)'(y-y)+(,L3-;3)'X'X(f3-f3)J} (8.5)

where ;3 = (X'X)- 1 X'y, andy= X)3. It then follows that ;3 is sufficient for ,13 if a- 2 is known, and (;3, (y-y)'(y-y)) is jointly sufficient for (,13, a- 2 ). Further,

and is independent of

We take the prior,

(8.6) This leads to the posterior,

242

8 Some Common Problems in Inference (8.7)

It can be seen that ,B[,B,a- 2

rv

Nv(,B,a- 2 (X'X)- 1 )

and that the posterior distribution of a- 2 is proportional to an inverse x;,-p" Integrating out a- 2 from this joint posterior density yields the multivariate t marginal posterior density for ,8, i.e., n(,B[y)

=

F(n/2)[X'X[ 112 s-P (r(1/2))P F((n- p)j2)(y'n- p)P

[

(,B-,B)'X'X(,B-,8)]-~

1+-'-----:-'---,---'::-----'-

(n- p)s2

'

(8.8)

where s 2 = (y- y)'(y- y)j(n- p). From this, it can be deduced that the posterior mean of ,B is ,B if n 2: p + 2, and the 100(1-a )% HPD credible region for ,B is given by the ellipsoid (8.9) where Fp,n-p(a) is the (1 -a) quantile of the Fp,n-p distribution. Further, if one is interested in a particular /3j, the fact that the marginal posterior distribution of /3j is given by (8.10) where djj is the jth diagonal entry of (X' x)- 1 , can be used. Conjugate priors for the normal regression model are of interest especially if hierarchical prior modeling is desired. This discussion, however, will be deferred to the following chapters where hierarchical Bayesian analysis is discussed. Example 8.1. Table 8.1 shows the maximum January temperatures (in degrees Fahrenheit), from 1931 to 1960, for 62 cities in the U.S., along with their latitude (degrees), longitude (degrees) and altitude (feet). (See Mosteller and Tukey, 1977.) It is of interest to relate the information supplied by the geographical coordinates to the maximum January temperatures. The following summary measures are obtained.

X' X

=

62.0 2365.0 5674.0 56012.0 2365.0 92955.0 217285.0 2244586.0 5674.0 217285.0 538752.0 5685654.0 [ 56012.0 2244586.0 5685654.0 1.7720873 X 108

l '

8.2 Linear Regression

243

Table 8.1. Maximum January Temperatures for U.S. Cities, with Latitude, Longitude, and Altitude City Latitude Longitude Altitude Max. Jan. Temp Mobile, Ala. 30 88 5 61 32 86 160 59 Montgomery, Ala. 58 134 50 30 Juneau, Alaska 33 112 1090 64 Phoenix, Ariz. 286 Little Rock, Ark. 34 51 92 340 Los Angeles, Calif. 34 118 65 San Francisco, Calif. 65 55 37 122 5280 Denver, Col. 42 39 104 40 New Haven, Conn. 41 72 37 Wilmington, Del. 41 135 39 75 44 Washington, D.C. 25 38 77 20 Jacksonville, Fla. 81 67 38 Key West, Fla. 74 24 5 81 Miami, Fla. 25 10 76 80 52 Atlanta, Ga. 84 1050 33 Honolulu, Hawaii 21 21 157 79 43 2704 Boise, Idaho 116 36 Chicago, Ill. 41 33 595 87 Indianapolis, Ind. 37 39 710 86 29 41 Des Moines, Iowa 805 93 27 42 620 Dubuque, Iowa 90 42 1290 Wichita, Kansas 97 37 Louisville, Ky. 44 450 38 85 64 New Orleans, La. 29 5 90 32 43 Portland, Maine 25 70 44 20 Baltimore, Md. 39 76 42 Boston, Mass. 21 37 71 42 Detroit, Mich. 585 83 33 84 23 Sault Sainte Marie, Mich. 46 650 Minneapolis -St. Paul, Minn. 815 22 44 93 40 St. Louis, Missouri 455 38 90 Helena, Montana 29 4155 112 46 Omaha, Neb. 32 1040 41 95 32 Concord, N.H. 290 43 71 Atlantic City, N.J. 43 74 10 39 Albuquerque, N.M. 46 106 4945 35 continues

(X' x)- 1 =

10- 5

94883.1914 -1342.5011 -485.0209 2.5756] -1342.5011 37.8582 -0.8276 -0.0286 -485.0209 -0.8276 5.8951 -0.0254 [ 2.5756 -0.0286 -0.0254 0.0009

244

8 Some Common Problems in Inference Table 8.1 continued

City Latitude Longitude Altitude Max. Jan. Temp 42 31 Albany, N.Y. 73 20 New York, N.Y. 40 40 73 55 Charlotte, N.C. 35 80 720 51 Raleigh, N.C. 35 78 365 52 Bismark, N.D. 46 1674 100 20 Cincinnati, Ohio 39 84 41 550 Cleveland, Ohio 41 81 35 660 Oklahoma City, Okla. 35 97 1195 46 45 Portland, Ore. 122 77 44 Harrisburg, Pa. 40 76 365 39 Philadelphia, Pa. 39 40 75 100 Charlestown, S.C. 32 61 79 9 Rapid City, S.D. 44 34 103 3230 Nashville, Tenn. 36 450 49 86 Amarillo, Tx. 35 101 3685 50 Galveston, Tx. 29 94 61 5 Houston, Tx. 29 64 95 40 Salt Lake City, Utah 40 111 4390 37 Burlington, Vt. 44 110 73 25 Norfolk, Va. 36 76 50 10 Seattle-Tacoma, Wash. 47 122 44 10 Spokane, Wash. 47 117 1890 31 Madison, Wise. 43 89 26 860 Milwaukee, Wise. 43 87 28 635 Cheyenne, Wyoming 41 104 6100 37 San Juan, Puerto Rico 18 66 35 81

I

X y =

2739.0 ) ( 100.8260) 99168.0 -1.9315 252007.0 'f3 = 0.2033 ' and s = 6.05185. A

(

2158463.0

-0.0017

On the basis of the analysis explained above, /3 may be taken as the estimate of (3 and the HPD credible region for it can be derived using (8.9). Suppose instead one is interested in the impact of latitude on maximum January temperatures. Then the 95% HPD region for the corresponding regression coefficient (32 can be obtained using (8.10). This yields the t-interval, ;32 ± sy'd22t 58 (.975), or (-2.1623, -1. 7007), indicating an expected general drop in maximum temperatures as one moves away from the Equator. If the joint impact of latitude and altitude is also of interest, then one would look at the HPD credible set for ((32 , (34 ). This is given by { (/32, /34) : (/32

+ 1.9315,/34 + 0.0017)C- 1 (/32 + 1.9315,/34 + 0.0017)' :'::: 2s 2F2,5s(a) },

8.3 Logit Model, Probit Model, and Logistic Regression

245

oM 0 0 0

oM 0 0 0 I

N 0 0 0 I

1"1 0 0 0 I

'
-2.4

-2

-2.2

-1.8

-1.6

-1.4

Fig. 8.1. Plot of 95% and 99% HPD credible regions for

where cis the appropriate 2 -4 C = 10

(

X

(/32 , /34 ).

2 block from (X' x)-1,

3.7858 -2.8636 X 10- 3 ) -2.8636 X 10- 3 9.2635 X 10- 5

.

Plotted in Figure 8.1 are the 95% and 99% HPD credible regions for

((32 , (34 ). Impact of altitude on maximum temperatures seems to be very limited for the case under consideration. Literature on Bayesian approach to linear regression is very large. Some of this material relevant to the discussion given above may be found in Box and Tiao (1973), Leamer (1978), and Gelman et al. (1995).

8.3 Logit Model, Probit Model, and Logistic Regression We consider a problem here that is related to linear regression but actually belongs to a broad class of generalized linear models. This model is useful for problems involving toxicity tests and bioassay experiments. In such experiments, usually various dose levels of drugs are administered to batches of animals. Most often the responses are dichotomous because what is observed is whether the subject is dead or whether a tumor has appeared. This leads to a setup that can be easily understood in the context of the following example.

246

8 Some Common Problems in Inference

Example 8.2. Suppose that k independent random variables y 1 , Y2, ... , Yk are observed, where Yi has the B(ni, Pi) probability distribution, 1 ::::; i ::::; k. Yi may be the number of laboratory animals cured of an ailment in an experiment involving ni such animals. It is certainly possible to make inference on each Pi separately based on the observed Yi (as discussed previously). This, however, is not really useful if we want to predict the results of a similar experiment in future. Suppose that the Pi are related to a covariate or an explanatory variable, such as dosage level in a clinical experiment. Then the natural approach is regression as described in the previous problem, because this allows us to explore and present the relationship between design (explanatory) variables and response variables, and (if needed) predictions of response at desired levels of the explanatory variables. Let ti be the value of the covariate that corresponds with Pi, i = 1, 2, ... , k. Because Pi's are probabilities, linking them to the corresponding ti 's through a linear map as was done earlier does not seem appropriate now. Instead it can be made through a link function H such that Pi = H((30 + {3 1 ti)· H, here, is a known cumulative distribution function (c.d.f.) and (30 and (31 are two unknown parameters. (If H is an invertible function, this is precisely H- 1 (pi) = (30 + {3 1 ti-) If the standard normal c.d.f. is used for H, the model is called the probit model, and if the logistic c.d.f. (i.e., H(z) = e-z /(1+e-z)) is used, it is called the logit model. The likelihood function for the unknown parameters, (30 and (3 1 , is then given by

Suppose 1r((30 , {3 1 ) is the prior density on ((30 , {3 1 ) so that the posterior density is

It may be noted that the sample size ni and the covariates ti (dose level) are

treated here as fixed, or equivalently the model conditional on those variables is analyzed (as in the linear regression problem). More generally, instead of a single covariate t, if we have a set of s covariates represented by x and the corresponding vector of coefficients {3, the posterior density of f3 will be k

1r(/3jdata) ex 1r(f3)

IT H(f3'xi)Y' (1- H(f3'xi)t'-y'.

(8.11)

i=1

8.3.1 The Logit Model

If we use the logit model whereby Pi= exp{ -f3'xi}/{1 + exp{ -f3'xi} }, and hence -log(pj(1- Pi))= f3'xi, we obtain the likelihood function

8.3 Logit Model, Probit Model, and Logistic Regression

247

which is largely intractable (but see Problem 7). Usually a flat prior such as

1r({3) ex: 1 is employed, but because f3 can be considered to be regression coefficients under reparameterization, an approximate conjugate normal prior can also be used. In such a case, a hierarchical prior structure is also meaningful. To motivate the approximate conjugate hierarchical prior structure, consider the following large sample approximation. For simplicity, let there be only one covariate t. Assume that the ni are large enough for a satisfactory normal approximation of the binomial model. If we let Pi = yi/ ni, then (approximately) these Pi are independent N (Pi, Pi ( 1 - Pi) I ni) random variates. Now let ei = -log(pi/(1- Pi)) and {Ji = -log(pi/(1- Pi)). It can be shown that, approximately, (ei - ei)Jnipi(1- Pi) are independent N(O, 1) random variates. Then (again approximately), the likelihood function for ((30 , (3 1 ) is (8.13) where Wi = niPi(1- Pi) are known weights. This suggests a bivariate normal prior on ((30 , (3 1 ) as the first level in the hierarchical structure. Now the problem is very similar to that discussed in Section 5.4. Further, there is also another important use for the approximate likelihood in (8.13). Its product with the conjugate normal prior discussed above can be used as a natural proposal density for the M-H algorithm in the computations required for inferential purposes (see Problem 9). If instead a flat prior on f3 is to be employed, then (8.13) itself (up to a constant) can be used as the proposal density.

Example 8.3. (Example 8.2 continued). The data given in Table 8.2 is from Finney (1971) (which originally appeared in Martin (1942)) where results of spraying rotenone of different concentrations on the chrysanthemum aphids in batches of about fifty are presented. The concentration is in milligrams per liter of the solution and the dosage x is measured on the logarithmic scale. The median lethal dose LD50, the dose at which 50% of the subjects will expire, is one of the quantities of inferential interest. The plot of p = yIn against x shown in Figure 8. 2 is S-shaped, so a linear fit for p based on x is unsatisfactory. Instead, as suggested by Figure 8.3, the logistic regression is more appropriate here. Suppose that a flat prior on ((30 , (31 ) is to be used. Then the implementation of M-H algorithm as explained above using the bivariate normal proposal density is straightforward. A scatter plot of 1000 values of ((30 , (3 1 ) so obtained is shown in Figure 8.4.

248

8 Some Common Problems in Inference

N 0

0.4

0.6

0.8

1.2

1

Fig. 8.2. Plot of proportion of deaths against dosage level.

N 0

0.2

0.4

0.6

0.8

Fig. 8.3. Plot of logistic response function:

1

1.2

e-5+?x /(1

+ e-5+?x).

8.3 Logit Model, Probit Model, and Logistic Regression

249

Table 8.2. Toxicity of Rotenone Concentration (Ci) Dose (xi= log(ci)) Batch Size (ni) Deaths (yi) 2.6 0.4150 6 50 16 48 3.8 0.5797 24 46 5.1 0.7076 42 49 7.7 0.8865 44 1.0086 50 10.2

"7 I 0

"""" \D I

-rl

I

fll

+' (!)

.Q

<:0 I

~

-rl I

3

4

5 beta 0

5

7

Fig. 8.4. Scatter plot of 1000 (!30, !31) values sampled using M-H algorithm.

A comparison of the estimates of {30 and {31 obtained using MLE and posterior means are shown in Table 8.3. Histograms of the M-H samples of {30 and /31 are shown in Figure 8.5 and Figure 8.6, respectively. They seem to be skewed and hence the posterior estimates seem more appropriate. Let us consider the estimation of LD50 next. Note that for the logit model LD50 is that dosage level t 0 for which E(yi/nilti = t 0 ) = exp( -({30 + {31 to))/(1 + exp(-(/30 + {31 to))) = 0.5. This means that LD50 = -/30 /{3 1 . We Table 8.3. Estimates of f)0 and !31 from Different Methods Method tJo s.e. (tJo) !31 s.e. (!31) MLE from logistic regression 4.826 -7.065 Posterior estimates 4.9727 0.6312 -7.266 0.8859

250

8 Some Common Problems in Inference

-

r--~

-

-

-

-

-

r-

r--

-

I--

-

-

~

f---1

Fig. 8.5. Histogram of 1000 f3o values sampled using M-H algorithm.

-

r-

1-1--

r-

r-

,---

r-r-

1--

-

-~-----

r-

1--

,--r-

~

-10

~

__r-8

-6

Fig. 8.6. Histogram of 1000 /3r values sampled using M-H algorithm.

-4

8.3 Logit Model, Probit Model, and Logistic Regression

251

can easily estimate this using our M-H sample, and we obtain an estimate of 0.6843 (in the logarithmic scale) with a standard error of 0.022.

8.3.2 The Probit Model As mentioned previously, if the standard normal c.d.f., is used for the link function H above, we obtain the probit model. Then, assuming that n(/3) is the prior density on f3 the posterior density is obtained as k

n(f31data) ex n(/3)

II <J>(f3'xi)Y; (1- <J>(f3'xi)t;-y;.

(8.14)

i=l Analytically, this is even less tractable than the posterior density for the logit model. However, the following computational scheme developed by Albert and Chib (1993) based on the Gibbs sampler provides a convenient strategy. First consider the simpler case involving Bernoulli y/s, i.e., ni = 1. Then, k

n(f31data) ex n(/3)

II <J>(f3'xi)Y; (1- <J>(f3'xi)) -y;. 1

i=l The computational scheme then proceeds by introducing k independent latent variables Z 1 , Z 2 , ... , Zk, where Zi '"'"' N(f3'xi, 1). If we let Yi = I(Zi > 0), then Y1 , ... , Yk are independent Bernoulli with Pi = P(Yi = 1) = <J>(f3'xi)Now note that the posterior density of f3 and Z = (Z1 , ... , Zk) given y = (Yl, · .. , Yk) is k

n(/3, Zly) ex n(/3)

II {I(Zi > O)I(yi = 1) + I(Zi:::; O)J(yi = 0)} ¢(Zi- f3'xi), i=l

(8.15) where¢ is the standard normal p.d.f. Even though (8.15) is not a joint density which allows sampling from it directly, it allows Gibbs sampler to handle it since n(f31Z, y) and n(ZI/3, y) allow sampling from them. It is clear that k

n(/3IZ,y) ex n(/3)

II ¢(Zi- f3'xi),

(8.16)

i=l whereas

N(f3'xi, 1) truncated at the left by 0 if Yi = 1; zi I(3, y'"'"' { N(f3'xi, 1) truncated at the right by 0 if Yi = 0.

(8.17)

Sampling Z from (8.17) is straightforward. On the other hand, (8.16) is simply the posterior density for the regression parameters in the normal linear model with error variance 1. Therefore, if a flat noninformative prior on f3 is used, then

252

8 Some Common Problems in Inference Table 8.4. Lethality of Chloracetic Acid

Dose Fatalities .0794 1 .1000 2 .1259 1 .1413 0 .1500 1 .1558 2

Dose Fatalities .1778 4 .1995 6 .2239 4 .2512 5 .2818 5 .3162 8

,BjZ,y"' Ns(/3z, (X'x)- 1 ), where X = (x~, ... , x~)' and /3z = (X'X)- 1 X'Z. If a proper normal prior is used a different normal posterior will emerge. In either case, it is easy to sample ,8 from this posterior conditional on Z. Extension of this scheme to binomial counts Y1 , Y2 , ... , Yk is straightforward. We let Yi = I:j~ 1 Yi 1 where Yi 1 = I(Z; 1 > 0), with Z;1 "'N(,B'x;, 1) are independent, j = 1, 2, ... , ni, i = 1, 2, ... , k. We then proceed exactly as above but at each Gibbs step we will need to generate I:i= 1 n; many (truncated) normals Z; 1 . This procedure is employed in the following example. Example 8.4. (Example 8.2 continued). Consider the data given in Table 8.4, taken from Woodward et al. (1941) where several data sets on toxicity of certain acids were reported. This particular data set examines the relationship between exposure to chloracetic acid and the death of mice. At each of the dose levels (measured in grams per kilogram of body weight), ten mice were exposed. The median lethal dose LD50 is again one of the quantities of inferential interest. The Gibbs sampler as explained above is employed to generate a sample from the posterior distribution of ((30 ,(31 ) given the data. A scatter plot of 1000 points so generated is shown in Figure 8. 7. From this sample we have the estimate of (-1.4426, 5.9224) for ((30 , (31 ) along with standard errors of 0.4451 and 2.0997, respectively. To estimate the LD50, note that for the probit model LD50 is the dosage level t 0 for which E(y;jn;jti = t 0 ) = (f3o+f31 t 0 ) = 0.5. As before, this implies that LD50 = - (30 / (3 1 . Using the sample provided by the Gibbs sampler we estimate this to be 0.248.

8.4 Exercises 1. Show how a random deviate from the Student's tis to be generated. 2. Construct the 95% HPD credible set for ()r -fh for the two-sample problem in Section 8.1 assuming = a~.

ai

8.4 Exercises

253

...-i

If)

I

-3

-2

-1

beta_O

Fig. 8. 7. Scatter plot of 1000 (f3o, /h) values sampled using Gibbs sampler.

3. Show that Student's t can be expressed as scale mixtures of normals. Using this fact, explain how the 95% HPD credible set for fh - (h can be constructed for the case given in (8.3). 4. Consider the data in Table 8.5 from a clinical trial conducted by Mr. S. Sahu, a medical student at Bangalore Medical College, Bangalore, India (personal communication). The objective of the study was to compare two treatments, surgical and non-surgical medical, for short-term management of benign prostatic hyperplasia (enlargement of prostate). The random observable of interest is the 'improvement score' recorded for each of the patients by the physician, which we assume to be normally distributed. There were 15 patients in the non-surgical group and 14 in the surgical group. Table 8.5. Improvement Scores

Apply the results from Problems 2 and 3 above to make inferences about the difference in mean improvement in this problem. 5. Consider the linear regression model (8.4). (a) Show that (!3, (y - y)' (y - y)) is jointly sufficient for ({3, 0" 2 ).

254

8 Some Common Problems in Inference

(b) Show that .i3ia 2 "'Nv(!3, a 2 (X' x)- 1 ), (y- y)'(y- y))ia 2 "'a 2 x;,_v, and they are independently distributed. (c) Under the default prior (8.6), show that .BIY has the multivariate t distribution having density (8.8). 6. Construct 95% HPD credible set for (/32 , {33 ) in Example 8.1. 7. Consider the model given in (8.12). (a) What is a sufficient statistic for ,8? (b) Show that the likelihood equations for deriving MLE of ,8 must satisfy

8. Justify the approximate likelihood given in (8.13). 9. Consider a multivariate normal prior on ,8 for Problem 7. (a) Explain how the M-H algorithm can be used for computing the posterior quantities. (b) Compare the above scheme with an importance sampling strategy where the importance function is proportional to the product of the normal prior and the approximate normal likelihood given in (8.13). 10. Apply the results from Problem 9 to Example 8.3 with some appropriate choice for the hyperparameters. 11. Justify (8.16) and (8.17), and explain how Gibbs sampler handles (8.15). 12. Analyze the problem in Example 8.4 with an additional data point having 9 fatalities at the dosage level of 0.3400. 13. Analyze the problem in Example 8.4 using logistic regression. Compare the results with those obtained in Section 8.3.2 using the probit model.

9

High-dimensional Problems

Rather than begin by defining what is meant by high-dimensional, we begin with a couple of examples.

=

Example 9.1. (Stein's example) Let N(J.Lpx 1 , Epxp a 2 Ipxp) be a p-variate normal population and X;= (X; 1 , ... , X;p), i = 1, ... , n ben i.i.d. p-variate samples. Because E = a 2 I, we may alternatively think of the data asp independent samples of size n from p univariate normal populations N(J-Lj, a 2 ), j = 1, ... ,p. The parameters of interest are the J-L/S. For convenience, we initially assume a 2 is known. Usually, the number of parameters, p, is large and the sample size n is small compared with p. These have been called problems with large p, small n. Note that n in Stein's example is the sample size, if we think of the data as a p-variate sample of size n. However, we could also think of the data as univariate samples of size n from each of p univariate populations. Then the total sample size would be np. The second interpretation leads to a class of similar examples. Note that the observations are not exchangeable except in subgroups, in this sense one may call them partially exchangeable. Example 9.2. Let f(x[J-Lj),j = 1, ... ,p, denote the densities for p populations, and X;j, i = 1, ... , n,j = 1, ... ,p denote p samples of size n from these p populations. As in Example 9.1, f(x[J-Lj) may contain additional common parameters. The object is to make inference about the /-Lj 's. In several path-breaking papers Stein (1955), James and Stein (1960), Stein (1981), Robbins (1955, 1964), Efron and Morris (1971, 1972, 1973a, 1975) have shown classical objective Bayes or classical frequentist methods, e.g., maximum likelihood estimates, will usually be inappropriate here. See also Kiefer and Wolfowitz (1956) for applications to examples like those of Neyman and Scott (1948). These approaches are discussed in Sections 9.1 through 9.4, with stress on the parametric empirical Bayes (PEB) approach of Efron and Morris, as extended in Morris (1983).

256

9 High-dimensional Problems

It turns out that exchangeability of J-Lr, ... , f.-Lp plays a fundamental role in all these approaches. Under this assumption, there is a simple and natural Bayesian solution of the problem based on a hierarchical prior and MCMC. Much of the popularity of Bayesian methods is due to the fact that many new examples of this kind could be treated in a unified way. Because of the fundamental role of exchangeability of J-lj 's and the simplicity, at least in principle, of the Bayesian approach, we begin with these two topics in Section 9.1. This leads in a natural way to the PEB approach in Sections 9.2 and 9.3 and the frequentist approach in Section 9.4. All the above sections deal with point or interval estimation. In Section 9.6 we deal with testing and multiple testing in high-dimensional problems, with an application to microarrays. High-dimensional inference is closely related to model selection in high-dimensional problems. A brief overview is presented in Sections 9.7 and 9.8. We discuss several general issues in Sections 9.5 and 9.9.

9.1 Exchangeability, Hierarchical Priors, Approximation to Posterior for Large p, and MCMC We follow the notations of Example 9.1 and Example 9.2, i.e., we consider p similar but not identical populations with densities f(xJJ-Lr), ... ,J(xiJ-Lp)· There is a sample of size n, X 1 j, ... , Xnj, from the jth population. These p populations may correspond with p adjacent small areas with unknown per capita income J-Lr, ... , f.-Lp, as in small area estimation, Ghosh and Meeden (1997, Chapters 4, 5). They could also correspond with p clinical trials in a particular hospital and J-Lj, j = 1, ... , p, are the mean effects of the drug being tested. In all these examples, the different studied populations are related to each other. In Morris (1983), which we closely follow in Section 9.2, the p populations correspond to p baseball players and J-Lj 's are average scores. Other such studies are listed in Morris and Christiansen (1996). In order to assign a prior distribution for the J-lj 's, we model them as exchangeable rather than i.i.d. or just independent. An exchangeable, dependent structure is consistent with the assumption that the studies are similar in a broad sense, so they share many common elements. On the other hand, independence may be unnecessarily restrictive and somewhat unintuitive in the sense that one would expect to have separate analysis for each sample if the J-lj 's were independent and hence unrelated. However, to justify exchangeability one would need a particular kind of dependence. For example, Morris (1983) points out that the baseball players in his study were all hitters; his analysis would have been hard to justify if he had considered both hitters and pitchers. Using a standard way of generating exchangeable random variables, we introduce a vector of hyperparameters 11 and assume J-L/S are i.i.d. 7r(J-LITJ) given '11· Typically, if f(xJJL) belongs to an exponential family, it is convenient to

9.1 Exchangeability, Hierarchical Priors, and MCMC

257

take n(t-ti1J) to be a conjugate prior. It can be shown that for even moderately large p-in the baseball example of Morris (1983), p = 18- there is a lot of information in the data on 1]. Hence the choice of a prior for 17 does not have much influence on inference about flj 's. It is customary to choose a uniform or one of the other objective priors (vide Chapter 5) for 1]. We illustrate these ideas by exploring in detail Example 9.1. Example 9.3. (Example 9.1, continued) Let f(xlt-tj) be the density of N(ftj, CT 2 ) where we initially assume CT 2 is known. We relax this assumption in Subsection 9.1.1. The prior for flj is taken to be N (77 1 , 772) where 77 1 is the prior guess about the flj 's and 772 is a measure of the prior uncertainty about this guess, vide Example 2.1, Chapter 2. The prior for 771,772 is 1r(77 1 ,772), which we specify a bit later. Because (Xj = 2:::~ 1 Xijln,j = 1, ... ,p) form a sufficient statistic and Xj 's are independent, the Bayes estimate for flj under squared error loss is E(t-t11X)

where X = (Xij, z Example 2.1)

= E(t-tjiX)

= 1, ... , n,

j

=

j

E(t-tjiX,1J)n(1JIX)d1J.

= 1, ... ,p),

X

= (X1, ... , Xp) and (vide

(9.1) with B = (CT 2In) I (772 + CT 2In), depends on data only through Xj. The Bayes estimate of flj, on the other hand, depends on Xj as above and also on (X 1 , . . . , Xp) because the posterior distribution of 17 depends on all the X1 's. Thus the Bayes estimate learns from the full sufficient statistic justifying simultaneous estimation of all the flj 's. This learning process is sometimes referred to as borrowing strength. If the modeling of flj 's is realistic, we would expect the Bayes estimates to perform better than the Xj 's. This is what is strikingly new in the case of large p, small n and follows from the exchangeability of flj 's. The posterior density n(1JIX) is also easy to calculate in principle. For known CT 2 , one can get it explicitly. Integrating out the flj 's and holding 1J fixed, we get Xj 's are independent and

(9.2) Let the prior density of (77 1 , 772) be constant on R x R +. (See Problem 1 and Gelman et al. (1995) to find out why some other choices like 1r(771,772) = 11772 are not suitable here.) Given (9.2) and n( 771, 772) as above,

258

9 High-dimensional Problems

-

1

"'p -

where X = P L..Jj=l XJ and S = In a similar way,

7r(JLJX)

=

"'p (XJ- - X)L..Jj=l

Jfi

2

.

7!'(/Lj IXj' 'TJ )7r( TJJX) dTJ,

(9.4)

j=l

where (!LJIXJ, TJ) are independent normal with .

.

mean as m (9.1) and vanance

'f/2a 'f/2

2

/n/

+a 2 n

(9.5)

and

(9.6) is the product of a normal and inverse-Gamma (as given in (9.3)). Construction of a credible interval for /Lj is, in principle, simple. Consider 7r(/LjJX) and fix 0