当前位置：文档库 › Bayesian Inference on Order-Constrained Parameters in Generalized Linear Models

Bayesian Inference on Order-Constrained Parameters in Generalized Linear Models

Bayesian Inference on Order Constrained Parameters in

Generalized Linear Models

David B.Dunson1,?and Brian Neelon2

July31,2002

1Biostatistics Branch,National Institute of Environmental Health Sciences,MD A3-03, P.O.Box12233,Research Triangle Park,NC27709

2Department of Biostatistics,University of North Carolina,Chapel Hill,NC

?email:dunson1@https://www.wendangku.net/doc/0a6782788.html,

S UMMARY.In biomedical studies,there is often interest in assessing the association be-tween one or more ordered categorical predictors and an outcome variable,adjusting for covariates.For a k-level predictor,one typically uses either a k?1degree of freedom(df) test or a single df trend test,which requires scores for the di?erent levels of the predic-tor.In the absence of knowledge of a parametric form for the response function,one can incorporate monotonicity constraints to improve the e?ciency of tests of association.This article proposes a general Bayesian approach for inference on order constrained parameters in generalized linear models.Instead of choosing a prior distribution with support on the constrained space,which can result in major computational di?culties,we propose to map draws from an unconstrained posterior density using an isotonic regression transformation. This approach allows?at regions over which increases in the level of a predictor have no ef-fect.Bayes factors for assessing ordered trends can be computed based on the output from a Gibbs sampling algorithm.Results from a simulation study are presented and the approach is applied to data from a time to pregnancy study.

K EY WORDS:Bayes factor;Constrained estimation;Categorical covariates;Gibbs sam-pling;Isotonic regression;Monotonicity;Simple ordering;Trend test.

1.Introduction

In analyzing data from epidemiologic studies investigating the relationship between a k-level ordered categorical predictor and an outcome variable,one typically uses a generalized linear model that includes indicator variables for the di?erent levels of the predictor and adjusts for covariates.Maximum likelihood estimates of the regression coe?cients can be obtained, and a k?1degree of freedom test can be used to assess overall evidence of an association between the predictor and the response.When following this strategy,analysts are often faced with noisy estimates of the category-speci?c regression coe?cients.This noise can lead to unreasonable patterns in the regression coe?cients corresponding to di?erent levels of the predictor,and to low power for the k?1degree of freedom test and for one degree of freedom contrasts.In order to improve e?ciency of overall tests for an association,one can take advantage of known orderings in the regression coe?cients.For example,in many applications,it is reasonable to assume that the regression coe?cients for the di?erent levels of a predictor follow a simple non-decreasing order.

This article proposes a Bayesian approach for incorporating order constraints on the regres-sion coe?cients in a generalized linear model.In the Bayesian paradigm,order constraints are typically imposed by choosing a prior distribution that has support on a restricted space. Compared with the vast literature on classical approaches for order restricted inference(cf., Robertson,Wright and Dykstra,1988),few Bayesian methods have been proposed.How-ever,as noted by Gelfand,Smith and Lee(1992),certain types of parameter constraints can easily be incorporated in Bayesian analyses through the use of truncated conjugate priors. In particular,it is straightforward to accommodate strict constraints(e.g.,θ1<θ2)in an analysis that uses Gibbs sampling for posterior computation by?rst choosing a prior den-sity without considering the constraint(e.g.,normal)and then discarding draws inconsistent with the constraint.Although such an approach is very useful in certain applications,it

does not allow for equalities in the parameters(e.g.,θ1=θ2),and hence is not appropriate for inferences on uncertain orderings.For example,one may be interested in assessing evi-dence of an increase in the probability of disease associated with increases in the level of an exposure.

To allow for?at regions over which increases in the level of a predictor have no e?ect on the response distribution,one could potentially choose a prior distribution for the category-speci?c regression coe?cients that allocates positive probabilities to equalities,such asθ1=θ2.However,such priors typically have non-conjugate structures and posterior computation may be di?cult,particularly in high dimensions.In fact,even in low dimensional problems, we have noted substantial di?culties in posterior computation.Markov chain Monte Carlo (MCMC)algorithms(cf.,Chen,Shao and Ibrahim,2000),such as the Metropolis-Hastings algorithm,tend to have high autocorrelation and low e?ciency due to the constraints on the parameters.Motivated by these di?culties and by frequentist approaches to order restricted inference,which apply an isotonic regression transformation to the unconstrained parameter estimates,this article proposes an alternative approach.

Following Gelfand et al.(1992),we?rst choose a prior for the regression coe?cients without considering the constraint.We then update the prior to obtain an unconstrained posterior density,which is available in closed form for conjugate models and can otherwise be estimated via Gibbs sampling(Gelfand and Smith,1990).We consider the constrained parameters as functions of the unconstrained parameters,calculated by applying the isotonic regression transformation.Under this framework,samples from the constrained posterior density can be obtained by transforming draws from the unconstrained posterior density.Hence,existing Gibbs sampling algorithms for posterior computation of GLMs can be used directly.

This approach is a Bayesian alternative to likelihood ratio(LR)tests constructed using

restricted maximum likelihood estimates(RMLEs)obtained under an order-constrained al-ternative hypothesis(cf.,Bartholomew,1961;Robertson et al.,1988;Agresti and Coull, 1996).As noted by Mukerjee and Tu(1993),the important problem of order-restricted inferences on regression coe?cients,which is our focus,has received relatively little atten-tion.Methods have been proposed for incorporating monotonicity constraints in the setting of generalized additive models(cf.,Bacchetti,1989;Morton-Jones et al.,2000),and these methods can be used for testing hypotheses about increasing response functions.Our goal is a general methodology for inferences on arbitrary order constraints on parameters in hi-erarchical models,with the immediate emphasis on regression parameters in GLMs.By following a Bayesian approach,it becomes straightforward to incorporate prior information and to perform inferences that do not rely on asymptotic approximations.

Section2proposes an approach for assessing evidence of ordering in a vector of normal means.Section3extends the approach to assess ordered trends in regression parameters in the setting of linear regression,probit models for categorical data,and generalized linear models.Section4applies the approach to data from a time to pregnancy study,and Section 5discusses the results.

2.Ordered Normal Means

2.1Isotonic Transformation Approach

Let y j denote an n j×1vector of iid N(μj,σ2)random variables,for j=1,...,k.Assuming a conjugate prior density,

θj|σ2～N(θ0j,σ2/κ0j),j=1,...,k,andσ2∝(σ2)?1,

the joint posterior density ofθ=(θ1,...,θk) andσ2is the product of

θj|σ2,y～N( θj, σ2j),j=1,...,k,andσ2|y～scaled Inv-χ2(νn,σ2n),(1)

where θj=(κ0jθ0j+n j y j)/(κ0j+n j), σ2j=σ2/(κ0j+n j),νn=n,and

νnσ2

n =

j=1

(n j?1)s2

κ0j n j

κ0j+n j

(y j?θ0j)2,

with s2

denoting the empirical variance of y j.

We are interested in an order constrained functional ofθ,which we denoteθ?=g(θ)∈Θ? k,where g(·)is the isotonic regression transformation.In particular,letting Σ=

diag( σ2

1,..., σ2

)denote the posterior covariance,we choose the following max-min formula

(Robertson et al.,1988;Hwang and Peddada,1994):

θ?j =g j(θ)=min

t∈U j

max

s∈L j

Σ?1

[s:t]

θ[s:t]

Σ?1

[s:t]

1t?s+1

,for j=1,...,k,(2)

with L j and U j denoting subsets of{1,...,k}such that the orderingθ?

j ≤θ?

is known for all

j ∈L j and the orderingθ?

j ≥θ?

is known for all j ∈U j.The[s:t]subscript indicates the

submatrices and subvectors corresponding to elements s,s+1,...,t.The isotonic regression transformation(2)setsθ?

equal to a weighted average of a subvector ofθ,where the weights depend on the posterior covariance Σand the elements of the subvector are chosen to satisfy the constraint.

The restricted maximum likelihood estimate(RMLE)is obtained by applying transformation (2)to the unrestricted maximum likelihood estimate,(y1,...,y k) ,withκ0j=0for j= 1,...,k(refer to Example1.2.1in Robertson et al.,1988).Formula(2)is a least squares projection from k to the restricted spaceΘ.We obtain draws from the posterior density for the restricted parameter,θ?,by transforming draws from the posterior density for the unrestricted parameter,θ,using a minimal distance mapping.In the limit as the sample size→∞,the posterior mode forθ?is equivalent to the RMLE.When the prior variance is large,we expect that the posterior mode will be similar to the RMLE even in small to moderate samples.

Suppose that interest focuses on assessing evidence to reject H0:θ?

1=θ?

=...=θ?

favor of H1:θ?

1≤θ?

≤...≤θ?

with at least one strict ordering.Note that H0and H1

are equivalent to H0:θ?

1=θ?

and H1:θ?

<θ?

under the constraint thatθ?∈Θ,with

Θ={θ?:θ?

1≤...≤θ?

}.As a measure of weight of evidence in favor of H1,we suggest the

Bayes factor(cf.,Kass and Raftery,1995),

BF10=Pr(H1|y)/Pr(H1)

Pr(H0|y)/Pr(H0)

Pr(θ?

<θ?

|y)/Pr(θ?

<θ?

)

Pr(θ?1=θ?

|y)/Pr(θ?1=θ?

)

,(3)

where Pr(θ?

1<θ?

)is the prior probability of H1and Pr(θ?

<θ?

|y)is the posterior prob-

ability.Unlike the posterior probability,the Bayes factor is adjusted for the prior odds of H1.Hence,BF10is a more direct measure of the weight of evidence in the current data,y, in favor of H1.The quantities needed to calculate the Bayes factor are derived in Appendix A.It is straightforward to show that BF10→∞for large samples and?nite prior precision under H1.

2.2Simulation Study

To study the behavior of the procedure,we conducted a small simulation study.In particular, we simulated data under3di?erent patterns inθ:

1.θ=(θ1,...,θk) =0[Null hypothesis,H0],

2.θ=(0,1/(k?1),2/(k?1),...,1) [Linearly increasing,H1],

3.θ=(0,...,0,1) [Threshold,H1].

For each pattern and for di?erent choices of n j(5,10,25),k(3,5)andκ0j(0.1,1.0),we simu-lated500data sets.Then,for each simulated data set,we calculated BF10assumingθ0=0. To assess the frequentist operating characteristics,we used BF10>0.95/0.05=19as the cuto?for signi?cance.This choice should produce results comparable to anα=0.05-level test.

The simulation results are reported in Table1,along with the results from a separate sim-

ulation evaluating the performance of the likelihood ratio(LR)test(cf.,Robertson et al., 1988,p63).Regardless of the values of k and n j and the choice ofκ0j,the type I error rates were close to the0.05nominal level for both the proposed Bayesian test and the LR test, with the error rates of the Bayes procedure slightly lower on average.The power of the two test procedures was similar overall,with the LR test having higher power under the linear alternative and the Bayes test having higher power under the threshold alternative.The power of the Bayes procedure was slightly higher when a more di?use prior was chosen(i.e., forκ0j=0.1).

3.Regression Models

3.1Linear Regression

It is straightforward to generalize the procedure described in Section2to perform inferences on orderings in regression parameters.Suppose data for each study subject consist of a response,y i,m ordered categorical predictors,c i1,...,c im,and an additional p predictors that are either continuous or binary,w i=(w i1,...,w ip) .For subject i(i=1,...,n),let

u il=(1(c

il=2),...,1(c

il=k l)

) denote a(k l?1)×1vector of dichotomous indicators for the level

of the l th categorical predictor,where c il∈{1,...,k l}.We focus on the linear regression model,

y i=β0+

l=1

βl+w

α+ i=x

θ+ i,(4)

whereβ0is an intercept parameter,βl=(βl1,...,βl,k

l?1

) are coe?cients for the di?erent

levels of c il,β=(β0,β

1,...,β

) ,αare coe?cients for w i,x i=(1,u

,...,u

,w i) ,

θ=(β ,α ) ,and i～N(0,τ?1).We assume conjugate priors forθandτ,θ～N(θ0,Σ0) andτ～gamma(a0,b0).

Let g Θl (·)denote the isotonic regression transformation mapping from

k l ?1→Θl ={β?l :0≤β?l 1≤...≤β?l,k l ?1},

and let β?l =g Θl (βl ).The transformation g Θl (·)produces a (k l ?1)×1vector,with elements de?ned in expression (2),but with a di?erent form for the posterior covariance (as described below).

Often,interest focuses on assessing evidence of an increasing trend in the mean response with increases in a categorical predictor,controlling for other factors.In particular,we may want to summarize evidence in the data for rejecting

H 0l :0=β?l 1=...=β?l,k l ?1(equivalent to β?l,k l ?1

=0for βl ∈Θl )in favor of the alternative hypothesis

H 1l :0≤β?l 1≤...≤β?l,k l ?1

with at least one strict inequality (equivalent to β?l,k l ?1

>0for βl ∈Θl ).Generalizing the procedure of Section 2,we recom-mend the Bayes factor,

BF l =Pr(H 1l |y ,X )/Pr(H 1l )Pr(H 0l |y ,X )/Pr(H 0l )=Pr(β?l,k l ?1>0|y ,X )/Pr(β?l,k l ?1>0)Pr(β?l,k l ?1=0|y ,X )/Pr(β?l,k l ?1=0),(5)

where X =(x 1,...,x n ) is the design matrix.

The Bayes factor,BF l ,can be estimated using the following procedure:

1.Draw samples from the posterior density of (θ,τ)using a Gibbs sampler,which alter-nately samples from the conditional densities,

[θ|τ,y ,X ]=N( θ, Σ)and [τ|θ,y ,X ]=gamma a 0+n 2,b 0+12n i =1(y i ?x i θ)2

,where θ= Σ(Σ?10θ0+τX y )and Σ=(Σ?10

+τX X )?1.

2.Draw samples from the prior density ofθby sampling from N(θ0,Σ0).

3.Apply the isotonic regression transformation to each sample from steps1and2to

obtain draws from the posterior and prior densities forβ?

,respectively(l=1,...,m).

4.Estimate the posterior and prior probabilities in expression(5)as proportions of the

draws forβ?

l,k l?1

equal to greater than0,for l=1,...,m.

3.2Probit Models for Categorical Data

Now suppose that y1,...,y n are independent Bernoulli distributed random variables,with

Pr(y i|x i,θ)=Φ

β0+

l=1

βl+w

=Φ(x

θ).(6)

As noted by Albert and Chib(1993),this is equivalent to assuming that y i=1(z

i>0)

,with z i～

N(x

θ,1)denoting independent and normally distributed random variables underlying y i,for i=1,...,n.Under the conditionally conjugate prior density forθde?ned in subsection3.1, posterior computation can proceed via a Gibbs sampler that alternates between(i)sampling from the conditional density of z i,

[z i|y i,x i,θ]=N(x

θ,1)truncated below(above)by0for y i=1(y i=0),

and(ii)sampling from the conditional density ofθ,[θ|z,y,X]=N( θ, Σ),where θand Σare as de?ned in subsection3.1,with z i substituted for y i andτ=1for identi?ability.The algorithm outlined in subsection3.1can be used to estimate order-constrained parameters and to compute Bayes factors for order restricted inferences.It is straightforward to extend this approach for categorical y i by following the approach of Albert and Chib(1993).

3.3Generalized Linear Models

Suppose that y1,...,y n are independent random variables from a distribution in the expo-nential family,withτdenoting a scalar dispersion parameter(possibly known)and with

canonical parameterψi related to covariates through the generalized linear model,

h(ψi)=ηi=β0+

l=1

βl+w

α=x

θ.(7)

Assuming a N(θ0,Σ0)prior density forθand a log-concave prior forτ,posterior computa-tion for(θ,τ)can proceed via an adaptive rejection Gibbs sampling algorithm(Dellaportas and Smith,1993).For normal linear regression and probit models,the unconstrained pos-terior covariance used to weight the unconstrained parameters in expression(2)is known. However,for other generalized linear models,such as logistic or log-linear regression models, the posterior covariance must be estimated.We recommend?rst drawing a large number of samples from the posterior density ofθand then setting Σequal to the empirical covari-ance.The estimate Σcan then be used within expression(2)to weight the elements ofθin mapping additional draws from the posterior density ofθontoΘ.

4.Application to Discrete Time to Event Data

4.1Time to Pregnancy Data and Model

We illustrate the methodology through application to time to pregnancy data from a study of dental assistants exposed to high levels of nitrous oxide(Rowland et al.,1992).Female dental assistants,aged19to39,were randomly selected from the dental-assistants registry of the California Department of Consumer A?airs and invited to participate in the study if they met eligibility criteria.These criteria included having a planned pregnancy within the past four years and having worked full time during the six months prior to having non-contracepting intercourse.Of the459eligible women who completed the screening questionnaire,428 provided detailed data on reproductive and contraceptive history,occupational exposures, and other factors related to fertility.Time to pregnancy was ascertained by calculating the number of menstrual cycles before the most recent pregnancy during which the woman was having non-contracepting sexual intercourse(cf.,Rowland et al.,1992).Summary statistics

are provided in Table2.

The goal of our analysis was to incorporate known order restrictions to improve e?ciency in assessing covariate e?ects on fecundability,the probability of conception in a menstrual cycle. We were particularly interested in the e?ect of smoking,which was borderline signi?cant[95% con?dence interval=(.27,1.06)]in a previous unconstrained maximum likelihood analysis (Rowland et al.,1992).To this end,we modelled the number of menstrual cycles until conception using a discrete-time hazard model that included23parameters measuring(i) the baseline hazard of conception for menstrual cycles1-13(individuals not conceiving by cycle13were censored);(ii)age;(iii)intercourse frequency;(iv)cigarettes smoked per day; and(v)use of oral contraceptives in the cycle prior to beginning the pregnancy attempt.

Let y=(y1,...,y n) denote the conception outcomes for the n=1968non-contracepting menstrual cycles under study,where y i=1if conception occurs in cycle i and y i=0 otherwise.We assume that

Pr(y i|x i,β)=Φ

β0+

l=1

βl

=Φ(x

θ),(8)

where u i1=(u i11,...,u i1,12) is a vector of0/1indicators of the cycle at risk for a women, u i2,u i3,u i4and u i5are vectors indicating,respectively,the category of age(19?24,25?29,> 30),intercourse frequency(<=1,1?3,3?4,>4times/week),cigarette smoking frequency (nonsmoker,1?5,6?10,11?15,>15cigs/day),and oral contraceptive use in the previous cycle(no,yes).

A Bayesian speci?cation of the model is completed with a conditionally conjugate N(θ0,Σ0) prior density for the regression parameters.The prior mean for the parameters,(β0,β1), characterizing the baseline time to pregnancy distribution was set equal to the posterior mean from an analysis of the Tietze(1968)data using model(8).The prior covariance for these parameters was set equal to the estimated posterior covariance from the Tietze

analysis,with the variance in?ated by a factor of10to accommodate possible di?erences between the study populations.The prior mean for the remaining parameters was set to0, while the prior covariance matrix was chosen to be diag(20,...,20).

We are interested in inferences on the order-constrained functional,θ?=gΘ(θ),of the unconstrained regression coe?cients,where

Θ={β?:0≥β?

11≥...≥β?

1,12

,0≤β?

≤β?

,0≥β?

≥...≥β?

,β?

≤0}.

The ordering inβ?

is consistent with published data from Tietze(1968),Wilcox et al.

(1988)and Rowland et al.(1992).Non-increasing elements ofβ?

imply that the discrete hazard of conception is non-increasing with increases in the number of menstrual cycles at risk.This re?ects the selection process by which women with relatively high fertility

conceive rapidly and are removed from later risk sets.The constraint onβ?

ensures that the hazard of conception is non-decreasing with increasing frequency of intercourse.Finally,the

constraints onβ?

4andβ?

restrict the hazard to be non-increasing with increasing number of

cigarettes smoked per day and with recent pill use.

4.2Results

Samples from the posterior density of the unconstrained regression coe?cients,θ,were ob-tained using the Gibbs sampler of Albert and Chib(1993),as described in subsection2.2. We ran the Gibbs sampler for5,000iterations,and discarded the?rst500as a burn-in. The isotonic regression transformation from 23→Θwas applied to each sample after the burn-in to obtain samples from the posterior density ofθ?.Posterior summaries of the order-constrained baseline conception rates are plotted in Figure1,along with posterior means for the unconstrained parameters and unconstrained maximum likelihood estimates. As expected,posterior means of the constrained parameters varied relatively smoothly with menstrual cycle number.There was no evidence of systematic di?erences between the es-

timates,suggesting that the order restriction is supported by the data.In addition,the posterior variance for the constrained parameters was41%lower(on average)than that for the unconstrained parameters,suggesting an improvement in e?ciency.

Table3presents posterior summaries of the regression coe?cients characterizing the co-variate e?ects.Table3also provides estimated Bayes factors for testing the hypotheses of interest;namely,decreases in the hazard of conception with increases in smoking frequency and with occurrence of recent pill use.For the highest category of cigarette smoking(>15 cigarettes smoked/day)there was strong evidence of decreased fecundability relative to the non-smokers,with the Bayes factor equal to28.76.This Bayes factor corresponds roughly to a p-value of0.03in terms of weight of evidence against the null hypothesis.In constrast, the p-value from the unconstrained maximum likelihood analysis presented by Rowland et al.(1992)was>0.05as was the Bayesian p-value(i.e.,posterior probability)from the un-constrained analysis.These results suggest an improvement in power for our order-restricted Bayesian procedure,since an association with smoking and prolonged time to pregnancy has been shown repeatedly in other studies.Both the constrained and unconstrained analyses showed clear evidence of lower fecundability for recent pill users,as anticipated.

To assess the sensitivity of the results to the choice of hyperparameters,we repeated the analyses for di?erent choices of prior.In particular,we considered priors forθwith(i) precision0.01instead of0.05;(ii)precision0.1instead of0.05;and(iii)with prior mean θ0=(0,...,0,?0.2,?0.5,0.2,0.4,0.6,?0.1,?0.2,?0.3,?0.4,?0.1) . Essentially identical results were obtained in each of our sensitivity analyses,with the pos-terior means(standard deviations)all within±0.02(0.01)of the values reported in Table 3.

5.Discussion

This article has proposed a general and easy-to-implement approach for assessing order-ings in regression parameters.Instead of de?ning a prior distribution with support on the constrained space,as in earlier methodology(cf.,Chen and Shao,1998;Gelfand and Kuo, 1991;Lavine and Mockus,1995;Ramgopal,Laud and Smith,1993),we consider the order-constrained parameters to be functions of parameters having support on k.A conjugate or conditionally-conjugate prior density is chosen and draws from the resulting unconstrained posterior density are mapped to the constrained space through use of an isotonic regression transformation.Posterior summaries of the order-constrained parameters and Bayes factors for testing of ordered trends can be calculated in closed form for the case of ordered normal means or directly from the output of Gibbs sampling algorithms for general regression prob-lems.Although we have focused on the problem of assessing simple ordering in regression coe?cients for di?erent categories of an ordinal predictor,the approach can be applied di-rectly for other types of orderings(e.g.,simple tree)and parameters(e.g.,scale parameters or random-e?ects).

For the simple problem of assessing ordering in normal means,the proposed Bayesian pro-cedure yields very similar inferences to a likelihood ratio test when a di?use prior is chosen. An advantage of the Bayesian procedure in this context is that exact posterior densities of the parameters can be easily estimated.In addition,it is straightforward to generalize the Bayesian procedure to arbitrary regression problems,while generalizations of the likelihood ratio test are limited by the di?culty of computing restricted maximum likelihood estimates in many cases.

Our approach is motivated primarily by convenience in computational implementation for a wide variety of problems.For example,with minimal programming e?ort in S-PLUS(Math-soft,Inc.,1999),Bayes factors and posterior summaries can be calculated directly from the

MCMC iterates output from WinBUGS(Lunn et al.,2000).Therefore,the procedure pro-vides a useful black box approach for Bayesian inference on order constrained parameters in hierarchical models.It should be noted that,from a Bayesian perspective,it is somewhat more natural to incorporate order constraints through an appropriately de?ned prior instead of through a transformation.However,posterior computation is often prohibitively di?cult when the prior density has constrained support,and general algorithms remain to be de-veloped.As we have demonstrated through a simulation study and through application to reproductive epidemiology data,our approach is a reasonable alternative.

A CKNOWLEDGEMENTS

The authors would like to thank Zhen Chen,Sean O’Brien and Shyamal Peddada for helpful comments and discussions.Thanks also to Andrew Rowland for providing the data for the example.

Appendix A

Derivation of Bayes Factor for Testing Ordering in Normal Means

From the de?nition of the isotonic regression transformation,it is apparent thatθ?

1=...=θ?

if and only ifθ1≥θ2≥...≥θk.Lettingδj=θj?θj+1,for j=1,...,k?1,the inequality δj≥0implies thatθj≥θj+1.The conditional posterior density ofδj|σ2is

δj|σ2,y～N( θj? θj+1, σ2j+ σ2j+1),j=1,...,k?1,

where σ2

=σ2/(κ0j+n j).Hence,conditional onσ2,the posterior probability ofδj≥0is

Pr(δj≥0|σ2,y)=Pr(θj≥θj+1|σ2,y)=Φ

θj? θj+1

σ2j+ σ2j+1

whereΦ(·)is the standard normal distribution function.Taking one minus the product across the di?erent elements of{δ1,...,δk?1},we obtain

1?Pr(θ1≥...≥θk|σ2,y)=Pr(H1|σ,y)=1?k?1

j=1

θj? θj+1

σ2j+ σ2j+1

A parallel approach can be used to derive the prior probability,Pr(H1|σ2).The Bayes factor can be calculated by plugging in the prior and posterior probabilities,conditional onσ2, into expression(3)and then removing the conditioning by numerical integration across the scaled-Inv-χ2(νn,σ2

)marginal posterior density forσ2.

References

Agresti,A.and Coull,B.A.(1996).Order-restricted tests for strati?ed comparisons of binomial proportions.Biometrics52,1103-1111.

Albert,J.H.and Chib,S.(1993).Bayesian analysis of binary and polychotomous response data.Journal of the American Statistical Association88,669-679.

Bartholomew,D.J.(1961).A test of homogeneity of means under restricted alternatives (with discussion).Journal of the Royal Statistical Society B23,239-281.

Chen,M.-H.and Shao,Q.-M.(1998)Monte Carlo methods on Bayesian analysis of con-strained parameter problems.Biometrika85,73-87.

Chen,M.-H.,Shao,Q.-M.and Ibrahim,J.G.(2000)Monte Carlo Methods in Bayesian Computation.New York:Springer.

Dellaportas,P.and Smith,A.F.M.(1993).Bayesian inference for generalized linear and proportional hazards models via Gibbs sampling.Applied Statistics42,443-459.

Gelfand,A.E.and Kuo,L.(1991)Nonparametric Bayesian bioassay including ordered poly-tomous response.Biometrika78,657-666.

Gelfand,A.E.and Smith,A.F.M.(1990)Sampling based approaches to calculating marginal densities.Journal of the American Statistical Association85,398-409.

Gelfand,A.E.,Smith,A.F.M.,and Lee,T.M.(1992)Bayesian analysis of constrained pa-rameter and truncated data problems using Gibbs sampling.Journal of the American Statistical Association87,523-532.

Hwang,J.T.G.and Peddada,S.D.(1994).Con?dence interval estimation subject to order restrictions.Annals of Statistics22,67-93.

Kass and Raftery(1995).Bayes factors.Journal of the American Statistical Association 90,773-795.

Lavine,M.and Mockus,A.(1995).A nonparametric Bayes method for isotonic regression.

Journal of Statistical Planning and Inference46,235-248.

Lunn,D.J.,Thomas,A.,Best,N.and Spiegelhalter,D.(2000).WINBUGS-A Bayesian modelling framework:Concepts,structure,and extensibility.Statistics and Computing 10,325-337.

Mathsoft,Inc.(1999).S-PLUS5for UNIX Guide to Statistics,Data Analysis Products Division,Mathsoft,Seattle.

McDermott,M.P.and Mudholkar,G.S.(1993).A simple approach to testing homogeneity of order-constrained means.Journal of the American Statistical Association88,1371-1379.

Morton-Jones,T.,Diggle,P.,Parker,L.,Dickinson,H.O.and Binks,K.(2000).Additive isotonic regression models in epidemiology.Statistics in Medicine19,849-859.

Mukerjee,H.and Tu,R.(1995).Order-restricted inferences in linear regression.Journal of the American Statistical Association90,717-728.

Ramgopal,P.,Laud,P.W.and Smith,A.F.M.(1993).Nonparametric Bayesian bioassay with prior constraints on the shape of the potency curve.Biometrika80,489-498. Robertson,T.,Wright,F.and Dykstra,R.(1988).Order Restricted Statistical Inference.

New York:Wiley.

Rowland,A.S.,Baird,D.D.,Weinberg,C.R.,Shore,D.L.,Shy,C.M.,and Wilcox,A.J.

(1992).Reduced fertility among women employed as dental assistants exposed to high levels of nitrous oxide.New England Journal of Medicine327,993-997.

Wilcox,A.J.,Weinberg,C.R.,O’Connor,J.F.,Baird,D.D.,Shlatterer,J.P.,Can?eld,R.E.

et al.(1988).Incidence of early loss of pregnancy.New England Journal of Medicine 319,189-194.

Tietze,C.(1968)Fertility after discontinuation of intrauterine and oral contraception.In-ternational Journal of Fertility13,385-389.

Results from simulation study of the frequentist operating characteristics of the proposed Bayes Factor approach for assessing ordering in normal means.

Rejection Rate

BF10Bayes?LR test?

Hypothesis Pattern k n jκ0j

H01350.1 1.570.070.07

H0135 1.0 1.540.05

H013100.1 1.440.040.05

H01310 1.0 1.460.04

H013250.1 1.600.030.05

H01325 1.0 1.490.03

H01550.1 1.810.050.04

H0155 1.0 1.680.02

H015100.1 1.770.030.04

H01510 1.0 1.800.03

H015250.1 1.930.040.06

H01525 1.0 1.710.01

H12350.117.20.480.42

H1235 1.010.90.40

H123100.152.90.680.65

H12310 1.046.60.65

H123250.1596.70.950.95

H12325 1.0506.90.94

H12550.114.10.440.44

H1255 1.08.550.28

H125100.130.20.600.77

H12510 1.021.50.54

H125250.1183.40.920.99

H12525 1.0135.50.90

H13350.131.20.590.43

H1335 1.016.50.48

H133100.1129.00.750.74

H13310 1.072.80.70

H133250.157740.980.98

H13325 1.067250.98

H13550.119.80.530.46

H1355 1.012.80.40

H135100.166.40.720.74

H13510 1.051.30.69

H135250.128570.990.99

H13525 1.020920.98

?proportion of simulations with BF10>19

?proportion of simulations with LR test p-value<0.05