On the frequentist and Bayesian approaches to hypothesis testing

Statistics & Operations Research Transactions SORT 30 (1) January-June 2006, 3-28 Statistics & Operations Research Transactions On the frequentist and Bayesian approaches to hypothesis testing c Institut d’Estad´ ıstica de Catalunya [email protected] ISSN: 1696-2281 www.idescat.net/sort El´ ıas Moreno and F. Javier Gir´ on University of Granada and University of M´alaga Abstract Hypothesis testing is a model selection problem for which the solution proposed by the two main statistical streams of thought, frequentists and Bayesians, substantially differ. One may think that this fact might be due to the prior chosen in the Bayesian analysis and that a convenient prior selection may reconcile both approaches. However, the Bayesian robustness viewpoint has shown that, in general, this is not so and hence a profound disagreement between both approaches exists. In this paper we briefly revise the basic aspects of hypothesis testing for both the frequentist and Bayesian procedures and discuss the variable selection problem in normal linear regression for which the discrepancies are more apparent. Illustrations on simulated and real data are given. MSC: 62A01; 62F03; 62F15 Keywords: Bayes factor, consistency, intrinsic priors, loss function, model posterior probability, p - values. 1 Introduction In parametric statistical inference estimating parameters and hypothesis testing are two basic but different problems, although some Bayesian estimation devices have gained some popularity as hypothesis testing tools. For instance, high posterior probability regions have been used as acceptance regions for the null hypothesis (see, for instance, the analysis of variance solution proposed by Lindley 1970, Box and Tiao 1992). This is misleading as far as the hypotheses to be tested do not play any role in the construction of such acceptance regions. Address for correspondence: El´ ıas Moreno and F. Javier Gir´ on. Department of Statistics and O.R. University of Granada and University of M´ alaga. Received: November 2005 4 On the frequentist and Bayesian approaches to hypothesis testing To distinguish between estimation and testing was strongly recommended by Jeffreys (1961, pp. 245-249), mainly because the methods commonly used for estimation are typically not suitable for hypothesis testing. Thus, we think it is timely to devote a separate paper to discuss the setting and tools devised for hypothesis testing from the frequentist and Bayesian perspectives. We want to thank the Editor of SORT for inviting us to contribute on this topic. In this paper we deal with frequentist and Bayesian parametric hypothesis testing procedures. A third approach, based solely on the likelihood function, which we do not discuss here, is to be found in Royall (1997) and Pawitan (2001). As Pawitan (2001, p. 15) states: The distinguishing view is that inference is possible directly from the likelihood function; this is neither Bayesian nor frequentist, and in fact both schools would reject such a view as they allow only probability-based inference. From the frequentist viewpoint two closely related methods have been developed. One is the Neyman-Pearson theory of significance tests and the other one is based on Fisher’s notion of p-values. Here, we shall give arguments that make the p-values to be preferable to significance tests. On the Bayesian side, robustness with respect to the prior showed that there is a strong discrepancy between the frequentist and Bayesian solutions to parametric hypothesis testing (Berger 1985, 1994, Berger and Delampady 1987, Berger and Sellke 1987, Berger and Mortera 1999, Casella and Berger 1997, Moreno and Cano 1998, Moreno 2005, among others). This means that the discrepancy is not due to the prior chosen for the Bayesian analysis but it is of a more fundamental nature which is inherent to the procedures. In particular, there is a marked difference on the way frequentist and Bayesian methods account for the sample size, and the dimensions of the null and the alternative parametric spaces. Since subjective prior distributions for the parameters of the models involved in hypothesis testing are not generally available, and their use is perceived as the weak point in the Bayesian implementation, objective prior distributions will be employed in this paper. By objective priors we mean priors that only depend on the sampling model and theoretical training samples. This priors are called intrinsic priors (Berger and Pericchi 1996, Moreno 1997, Moreno et al. 1998) and their merits can be judged for each specific application. We remark that they have been proved to behave extremely well in a wide variety of problems (Casella and Moreno 2004, 2005, Gir´ on and Moreno 2004, 2005, Moreno et al. 1999, 2000, 2003, 2005, Moreno and Liseo 2003). The rest of the paper is organized as follows. The second section is devoted to briefly describing significance tests and p-values. Section 3 reviews the Bayesian testing machinery and justifies the need for objective priors. The challenging problem of testing whether the means of two normal distributions with unknown variances are equal is considered in Section 4. A comparison of the frequentist and Bayesian testing procedures for the normal linear regression model, including a discussion on some fundamental issues of the variable selection problem, is given in Section 5, and some conclusions and recommendations are given in Section 6. El´ ıas Moreno, F. Javier Gir´ on 5 2 Significance tests and p-values Let Xdenote an observable random variable and x=(x1,...,xn)anavailable sample from either the model P0(x)orP1(x) with probability densities f0(x)orf1(x), respectively. Suppose that we want to choose between either the null hypothesis H0: f0(x) or the alternative H1:f1(x). This is the simplest hypothesis testing formulation. For this problem the well-known Neyman-Pearson theory of significance tests proposes a subset of the sample space Rn, the so-called critical or rejection region, Wα=y:f1(y) f0(y)≥cα,(1) as the region containing evidence against the null hypothesis. The threshold cαis determined so that the probability of the critical region under the null is α,thatis P0(Wα)=Wα f0(y)dy=α. Given the data x, the null H0is rejected at the significance level αif x∈Wα,and accepted otherwise. The value αis usually taken to be small so that we have a small probability of rejecting the null when it is true. Typically αis chosen to be 0.05 or 0.01. Fisher criticized the notion of significance of the Neyman-Pearson theory and proposed replacing cαin (1) with the observed likelihood ratio λn(x)=f1(x)/f0(x), so that the above critical region becomes Wn(x)=y:f1(y) f0(y)≥λn(x), that now depends on the observed data x. The probability of this region under the null, p=P0(Wn(x)) say, is called the p-value of the data x, and it is interpreted in such a way that a small enough p-value implies evidence against H0. The smaller the p-value the stronger the evidence against the null. Significance tests and p-values substantially differ. While the p-value depends on the data thus giving a measure of how strongly the observed data reject the null, the significance test does not provide such a measure. This simple fact implies that, for the dicotomous testing problem, the p-value is consistent under the null and the alternative, while the significance test is consistent under the alternative but it is not under the null. Indeed, when sampling from the null and the sample size grows to infinity the likelihood ratio statistic λn(x) tends to zero under P0, and the corresponding sequence of p-values tends to 1. This implies that, asymptotically, there is no evidence at all to reject the null. Alternatively, when the sample comes from the alternative hypothesis, the statistic λn(x) tends to ∞under P0, and hence the corresponding sequence of p- 6 On the frequentist and Bayesian approaches to hypothesis testing values tends to 0 showing an increasing evidence to reject the null. On the other hand, when sampling from the null, the significance test will reject H0with probability αeven when the sample size increases to infinity. This means that the significance test is not a consistent procedure under the null, although it is consistent when sampling from the alternative (Casella and Berger 1990, pp. 381-382, Wilks 1962, pp. 411). Under some regularity conditions on the likelihood, this theory can be extended to the more general situation of considering a family of probability densities {f(x|θ),θ ∈ Θ}, where Θmight be a multidimensional parametric space, and the null H0:θ∈Θ0 and the alternative H1:θ∈Θ1may contain more than one density. To preserve the good properties of the mentioned simplest hypothesis testing procedures, the hypotheses must be nested, that is Θ0⊂Θ1, and also we must have k0=dim(Θ0) strictly smaller than k1=dim(Θ1). The critical region now is Wn(x)=y:f(y|ˆ θ1(y)) f(y|ˆ θ0(y)) ≥ˆ λn(x), where ˆ θ1(y), ˆ θ0(y) are the MLE’s of θin the spaces Θ1and Θ0, respectively, and ˆ λn(x)=f(x|ˆ θ1(x))/f(x|ˆ θ0(x)). Two important difficulties now arise with the p-values. First, it is not clear how the probability of Wn(x) under the null can be computed when either the null is composite or the distribution induced by ˆ λn(y) depends on some (nuisance) parameter, as in the Behrens-Fisher problem which we briefly deal with in Section 4. Second, while a small p-value might contain evidence against the null for a dataset, the same p-value for a different sample size and/or for a null parameter space of different dimension might not contain the same evidence against the null. In fact, when we consider multiple tests, i.e. when we want to test that some subsets of regression coefficients are zero, as in variable selection, the frequentist literature (see, for instance Miller 2002) recognizes that the meaning of the p-value should depend on the dimensions of the null and the alternative parameter spaces and, consequently, provides a variety of methods to correct the p-value to account for this. It has also been recognized that the bigger the sample size the smaller the evidence of a given p-value against the null so that it is also necessary to correct the p-value to account for the sample size. These considerations have prompted the need to introduce frequentist criteria based on statistics that, in some way, adjust for all the varying parameters in the model such as the sample size and the dimensions of the null and the alternative hypothesis, such as the adjusted R2, Mallow‘s Cpor the AIC criteria. In summary, the meaning of the p-value is unfortunately unclear. Its interpretation should depend on the dimension of the null space, the dimension of the alternative space, and the sample size in an unknown, and probably complex and non-trivial, way; as a consequence, the calibration of a p-value is deemed to be a very difficult task, although some attempts for calibrating the p-values can be found in Sellke et al. (2001) and Gir´ on et al. (2004). El´ ıas Moreno, F. Javier Gir´ on 7 Furthermore, although a p-value is derived as a probabilistic statement and, consequently, lies between zero and one, it cannot be interpreted as the posterior probability that the null is true – this is an instance of the well known transposed conditional or prosecutor’s fallacy. However, practitioners of the p-values have very often this wrong probability interpretation in mind, maybe because this provides them with some sort of a (wrong) measurement devise for calibration. 3 Bayesian hypothesis testing From a Bayesian viewpoint the testing problem is treated as follows. Consider the simplest testing problem where we have to choose either the model M0:f0(x)or M1:f1(x) based on the observations x=(x1,...,xn). Let didenote the decision of choosing Miand let Pbe the prior probability defined on the model space {M0,M1}. Assume that a loss L(di,Mj)=cij,i,j=0,1, is incurred when we make the decision di and the true model is Mj(for other loss functions in model selection, see San Martini and Spezzaferri 1984, and Bernardo and Smith 1994). Assuming that the loss for a correct decision is zero, that is cii =0, and cij >0 otherwise, d1is the optimal decision when the posterior risks satisfy R(d0|x)>R(d1|x). This implies the following inequality P(M0|x)<c01 c01 +c10 . By Bayes theorem, this is equivalent to the inequality f1(x) f0(x)>c10 c01 P(M0) P(M1). Notice that the value of P(M0|x) is a measure, in a probability scale, of the strength we have in accepting the model M0. We now observe a first important difference between the p-value and the Bayesian report. While the former is obtained by integration over a region of the sample space, namely the rejection region, the latter is obtained directly from the loss function and the prior probabilities assigned to the models. The extension to the more realistic case of a parametric families of densities is straightforward. Consider the sampling models {f(x|θ0),θ 0∈Θ0}and {f(x|θ1),θ 1∈Θ1}. A complete Bayesian specification of the models needs prior distributions for the parameter θ0and θ1,thatis M0:{f(x|θ0),π 0(θ0)}, 8 On the frequentist and Bayesian approaches to hypothesis testing and M1:{f(x|θ1),π 1(θ1)}. Then under the loss function given above, assuming c01 =c10 and P(M0)=P(M1)= 1/2, the model M0is to be rejected if P(M0|x)=1 1+B10(x)<1/2,(2) where B10(x), the Bayes factor for models {M1,M0}, is the ratio of the marginal density, sometimes called the integrated or marginal likelihood, of the data under the two models, that is B10(x)=m1(x) m0(x)=f(x|θ1)π1(θ1)dθ1 f(x|θ0)π0(θ0)dθ0 . We note that a second important difference between the p-values and the Bayesian method is that while in the p-value approach the parameters of the null and the alternative hypothesis are estimated using the maximum likelihood method, in the Bayesian approach they are integrated out using prior distributions. For nested models it can be shown that, under mild conditions, the Bayesian procedure chooses the correct model with probability that tends to one as the sample size increases, so that it is a consistent procedure under both the null and the alternative hypothesis (see, for instance, O’Hagan and Forster 2004, p.182). This approach needs the specification of the losses cij,i,j=0,1, the prior distributions for parameters π0(θ0)andπ1(θ1), and the model prior (P(M0),P(M1)). While the specification of the loss function and the model prior seems to be a plausible task in real applications, the specification of subjective priors for parameters is admittedly a hard task. For instance, when θ1represents the vector of regression coefficients and the variance error of a regression model, to specify the prior is far from trivial. More so when the null parameter θ0is a subvector of the set of regression coefficients of θ1, thus indicating that a submodel is being considered plausible and a testing problem is called for. This problem is an important example in which the use of objective priors is fully justified. Unfortunately, the priors considered for estimation, as the Jeffreys or the reference priors by Berger and Bernardo (1992), are typically improper so that they depend on arbitrary multiplicative constants that leave the Bayes factor ill-defined as the following simple example shows. Example 1 Suppose that Xis a random variable with distribution N(x|μ, σ2), with both parameters unknown, and we want to test H0:μ=0 versus H1:μ0. This is equivalent to choosing between models M0:{x|σ0∼N(x|0,σ 0),π 0(σ0)}and El´ ıas Moreno, F. Javier Gir´ on 9 {x|μ, σ1∼N(x|μ, σ1),π 1(μ, σ1)}. The reference prior for the parameter of the null model is π0(σ0)=c0/σ0, where c0is an arbitrary positive constant that cannot be specified because π0is improper. Likewise, the reference prior for the parameter of the alternative model is π1(μ, σ1)=c1/σ1, where again c1is an arbitrary positive constant. Therefore, the Bayes factor B10(x)isdefined up to the multiplicative constant c1/c0, whatever the data xmight be. 3.1 Intrinsic priors Consider the Bayesian models M0:{f(x|θ0),π N 0(θ0)}, and M1:{f(x|θ1),π N 1(θ1)} where πN 0and πN 1are objective, or default, improper priors. Lempers (1971, section 5.3) overcomes the difficulty that the Bayes factor for improper priors is not well defined by considering a partial Bayes factor. This is a Bayes factor constructed as follows. A part of the sample x, the training part, is devoted to converting the reference improper priors for the parameters of the null and alternative models into proper posteriors. Then, with the rest of the sample the Bayes factor is computed using these proper posteriors as priors. That is, the whole sample xis split into (x=x(t),x(n−t)), where x(t) is the training sample of the vector and x(n−t)the remaining one. Then, the posterior distributions for the above priors are πi(θi|x(t)) =f(x(t)|θi)πN i(θi) f(x(t)|θi)πN i(θi)dθi ,i=0,1, that are now proper. Now, using the above posteriors as priors, the Bayes factor for the rest of the data x(n−t) turns out to be BP 10(x)=f(x|θ1)πN 1(θ1)dθ1 f(x|θ0)πN 0(θ0)dθ0f(x(t)|θ0)πN 0(θ0)dθ0 f(x(t)|θ1)πN 1(θ1)dθ1 =BN 10(x)B01(x(t)). Note that the partial Bayes factor is a well defined Bayes factor that uses each of the components of the sample only once. However, it does depend on the specific training sample x(t) we choose. 10 On the frequentist and Bayesian approaches to hypothesis testing To avoid the arbitrariness in choosing the training sample Berger and Pericchi (1996) suggested computing the partial Bayes factors for all possible training samples with minimal size x(), and then computing the mean of those partial Bayes factors. The number is chosen so that the marginals of x() with respect to the improper priors are positive and finite so that the factor B01(x()) is well-defined up to the multiplicative constant c0/c1. The resulting value was called the arithmetic intrinsic Bayes factor, which does not depend on any arbitrary constant nor the particular training sample. Then, the Bayes factor appearing in (2) is replaced with the arithmetic intrinsic Bayes factor for computing the null posterior probability “as if” it were a Bayes factor. We observe that the arithmetic intrinsic Bayes factor is a mean of (partial) Bayes factors and hence it reuses the sample observations. We also observe that for some subsamples of minimal size it might be the case that the marginal mi(x()) could be zero or infinite. In that case, the Bayes factor is not well defined and we adopt the convention of not considering those subsamples. This implies that the arithmetic intrinsic Bayes factor might be quite unstable depending on the nature of the sample at hand. Some samples can have very many nice subsamples of minimal size but others may not have so many. However, to use the arithmetic intrinsic Bayes factor “as if” it were a Bayes factor is, in our opinion, not the best use we can give to the arithmetic intrinsic Bayes factor. It can be better employed as a tool for constructing priors. In fact, the arithmetic intrinsic Bayes factor is not a Bayes factor although as the sample size increases it becomes more and more stable and tends to be a Bayes factor for the so called intrinsic priors. Thus, if we use theoretical training samples instead of actual samples along with a limiting procedure we end up with intrinsic priors (Moreno et al. 1998). Given the model M0:{f(x|θ0),π N 0(θ0)} and M1:{f(x|θ1),π N 1(θ1)}, where f(x|θ0) is nested into f(x|θ1)andπN 1is improper, the following statements can be proven. (i) The intrinsic prior for θ1conditional on an arbitrary but fixed point θ0is given by πI(θ1|θ0)=πN 1(θ1)EX()|θ1 f(X()|θ0) f(X()|θ1)πN 1(θ1)dθ1 , where X() is a vector of dimension with i.i.d components and distribution f(x|θ1), such that El´ ıas Moreno, F. Javier Gir´ on 11 0<f(X()|θ1)πN 1(θ1)dθ1<∞, being the smallest natural number satisfying the above inequality. Roughly speaking, coincides with the dimension of θ1. (ii) πI(θ1|θ0) is a probability density for θ1, for any fixed θ0. (iii) If the default prior πN 0(θ0) is also improper, the ratio BI 10(x)=f(x|θ1)πI(θ1|θ0)πN 0(θ0)dθ0dθ1 f(x|θ0)πN 0(θ0)dθ0 (3) is the limit of the sequence of Bayes factors given by Bi=Cif(x|θ1)πI(θ1|θ0)πN 0(θ0|Ci)dθ0dθ1 Cif(x|θ0)πN 0(θ0|Ci)dθ0 , where {Ci,i≥1}is a covering monotone increasing sequence of sets in Θ0.Of course it can be shown that the limiting value (3) does not depend on the chosen sequence {Ci,i≥1}. In summary, intrinsic priors are well defined priors for testing problem involving nested models. For some particular non-nested models intrinsic priors can also be defined (Cano et al. 2004). The Bayes factor for intrinsic priors can be seen as the stabilized version of the arithmetic intrinsic Bayes factor. Further, as the sample size n tends to infinity the sequence of intrinsic posterior probabilities of model M0 P(M0|x1,...,xn)=1 1+BI 10(x1,...,xn) tends to one when sampling from the null and tends to zero when sampling from the alternative, so that the intrinsic Bayesian procedure is consistent; for a result in this direction see Moreno and Gir´ on (2005a) for the case of the general normal linear model. 4 The two sample problem Ap-value does not always exist, so that some sort of “approximation” is in that case necessary. A classical example in which this situation ocurrs is that of comparing the means of two normal distributions with unknown variances. Let N(x1|μ1,σ 2 1), N(x2|μ2,σ 2 2) be two normal distributions where the means μ1,μ2and variances σ2 1,σ 2 2 are unknown. Suppose that samples x1=(x11,x12,...,x1n1)andx2=(x21,x22,...,x2n2) have been drawn independently for each of the distributions, and we are interested in testing the null H0:μ1=μ2versus the alternative H1:μ1μ2. 18 On the frequentist and Bayesian approaches to hypothesis testing where B10 =2(k+1)k0/2 ππ/2 0 sink0ϕ(n+(k+1) sin2ϕ)(n−k)/2 (nBn+(k+1) sin2ϕ)(n−k1)/2dϕ. (4) From this expression, and also from the frequentist analysis of the testing problem in subsection 5.1, it follows that for fixed values of the sample size n, the number of covariates kand the dimension of the null hypothesis k1, the statistic Bnis a sufficient statistic for the testing problem, as the Bayes factor for the intrinsic priors does not depend on other ancillary statistics such as happens with other Bayes factors for linear models found in the literature, which depend on the quotient of the determinants |XtX| and |Xt 1X1|. Bayesian testing procedures different from the above one have been given by Berger and Pericchi (1996), O’Hagan (1995) and Zellner (1986) who proposed the use of the arithmetic intrinsic Bayes factor, the fractional Bayes factor and the Bayes factor derived from the g-priors, respectively. Except for the arithmetic intrinsic Bayes factor, the other two proposals depend on some tuning parameters which have to be adjusted. Let us mention that for normal linear models the O‘Hagan fractional Bayes factor provides sensible fractional priors for testing problems in a similar asymptotic way as the arithmetic intrinsic Bayes factor provides intrinsic priors (Moreno 1997). For so doing, the tuning parameter in the fractional Bayes factor is fixed as the quotient m/n, where mis the minimal training sample size. The results obtained when using Bayes factors for fractional priors are very close to those provided by Bayes factors for intrinsic priors, and hence only intrinsic priors are being considered here. 5.3 Comparing the frequentist and Bayesian tests The dependence of the p-value and the posterior probability of the null on the sufficient statistic (Bn,n,k,k1), where n,k,andk1are ancillary makes possible the comparison of the frequentist and objective Bayesian test. For fixed values of the ancillaries n,kand k1,thep-value and the intrinsic posterior probability of the null model P(M0|Bn,n,k,k1) are monotone increasing functions of the Bnstatistic. This permits us to establish a one-to-one relation between both measures of evidence through the parametric equations y=P(M0|b,n,k,k1) p=Ibn−k 2,k−k1 2, (5) where the parameter b,thesufficient statistic, ranges in the interval [0,1]. El´ ıas Moreno, F. Javier Gir´ on 19 The separate behaviour of yand pas the sufficient statistic bgoes to zero or one is as follows. The null posterior probability and the p-value go to zero as btends to zero, whatever the values of the ancillaries, as lim b→0P(M0|b,n,k,k1)=0 and lim b→0Ibn−k 2,k−k1 2=0. If Bn=0, then the residual sum of squares of the full model SS is also 0; this means that there is no uncertainty in the full model, i.e. it is deterministic; thus, the reduced model M0has zero posterior probability and the p-value is also zero, so that the full model is obviously accepted. When Bntends to one, then the p-value tends to one whatever the values of the ancillaries, but the null posterior probability tends to a number strictly smaller that one (Theorem 2.2 in Gir´ on et al 2004) which, on the other hand, tends to one as ngoes to infinity. In this case, as Bn=1, the residual sum of squares of the full SS and the reduced null model SS1are the same so that the data favour either model equally; however, the frequentist evidence in favour of the reduced model M0is one as if there were no uncertainty about what model to choose but, on the other hand, the Bayesian test accounts for the uncertainty inherent in the data rendering a posterior probability of M0 greater than 1/2 but strictly less than 1; hence, the Bayesian test chooses the simpler model, which is a consequence of the built-in Occam’s razor implicit in the objective Bayesian test. From the above equations (5) we can eliminate the parameter bto obtain a explicit equation of the null posterior probability as a function of the p-value, n,kand k1, y=P(M0|I−1 p(n−k/2,k−k1/2),n,k,k1). From this equation it follows that for fixed values of n,kand k1, the null posterior probability is an increasing function of the p-value. Therefore, the difference between the frequentist and Bayesian measures of evidence is a simple calibration problem. For this reason this curve is given the name of calibration curve in Gir´ on et al (2004). Unfortunately, the null posterior probability also depends on n,k,k1so that the properties of the calibration curve have to be established in a case-by-case basis. When simultaneous hypothesis testing are considered we have to jump among different calibration curves hence losing the monotonicity between the null posterior probabilities and the p-values. In this latter case calibration is no longer possible. In the remainder of this section the number of possible regressors kwill be keep fixed. In Figures 1 and 2 we display the typical behavior of the calibration curves, first, for different values of the sample size nand fixed k1, and second, for different values of k1and fixed n. 20 On the frequentist and Bayesian approaches to hypothesis testing The calibration curves in Figure 1 correspond to k=10 and k1=9andn=20, 50 and 90; they indicate that for a given p-value the null posterior probability increases as the sample size increases. Further, from the consistency of the Bayes factor for intrinsic priors it follows that the slope of the calibration curve at the origin tends to infinity as nincreases. This shows that the evidence against the null conveyed by the p-values should be diminished as the sample size nincreases in order to reconcile the frequentist and Bayesian test. Otherwise, we would reject a null hypothesis that has a very high posterior probability. Figure 1:Three calibration curves for different sample sizes n =20,n=50 and n =90, when the number of regressors is k =10, and k1=9. Pr(H0|p) The curves in Figure 2, for a fixed value of the sample size n, indicate that small p-values correspond to small posterior probabilities when k1is large, and also that the posterior probability increases as k1decreases. This implies that for small values of k1a p-value would reject a null hypothesis that has a large posterior probability, a fact which is generally acknowledged in the literature. But, on the other hand, for large values of k1 ap-value would accept a null hypothesis that has a small posterior probability. Notice, in Figure 2, that for the curve with k1=9 there is an interval of p-values larger than 0.05 whose corresponding posterior probabilities of the null are smaller than 1/2. This important fact, which is generally overlooked by the feeling that p-values tend to reject the null hypothesis more often than the Bayesian tests, also reveals that the opposite may happen sometimes; namely, that a null hypothesis may be not rejected by the frequentist p-value and rejected by the Bayesian test for the same data. Figure 2 also indicates that, when comparing several nulls of different dimension k1that convey the same frequentist evidence, the Bayesian test chooses the simplest El´ ıas Moreno, F. Javier Gir´ on 21 Figure 2:Four calibration curves, plotted for values of p-values in the interval (0,0.10), for different choices of k1=0,5,7,9, when the sample size is n =50, and the number of regressors is k =10. Pr(H0|p) hypothesis or model. As before, this is another instance of the automatic Occam’s razor property implicit in the objective Bayesian test. 5.4 Variable selection in regression Although the purpose of this paper is to critically revise the similarities and, above all, the enormous differences between the two approaches to hypothesis testing, we want to devote this last section to the important problem of variable selection in normal linear regression, where we have to consider simultaneously a large number of multiple tests of different dimensions k1, and then provide an ordering of the plausibility of the different models considered. A recent contribution to the subject, mainly from a frequentist perspective, is the revised version of the classical monograph by Miller (2002), where a chapter on Bayesian variable selection has been added. We do not discuss here sequential, or non-exhaustive, frequentist variable selection procedures such as forward, backward and stepwise selection methods nor the more recent ones such as the lasso or shrinkage methods, in order to concentrate our discussion on exhaustive search methods from the frequentist and Bayesian perspectives. Neither do we discuss here the well known asymptotic BIC criterion for model selection, because we can dispense with it as we have the non-asymptotically based from above and from below Bayesian criteria. The all subsets selection criterion is based on the idea of classifying the set of all models, say M,inkdisjoint classes, where kis the number of covariates of the full model, according to the number of covariates j,i.e.M=k j=1Mj. Within each class, 22 On the frequentist and Bayesian approaches to hypothesis testing the model with minimum residual sum of squares SSjis chosen, so that we end up with k maximal models, one for each class. Note that p-values or the equivalent Fstatistic provide the same ordering within the class Mjand, consequently, the same sets of maximals, but these criteria turn to be useless when comparing models with different number of covariates, a fact which is well recognized by frequentist statisticians. Hence, the problem of choosing the best model within the maximals, is far from trivial and a large number of procedures have been proposed in the literature. The underlying idea is to correct the SSj, or a simple function of it and the ancillaries to account for the different number of covariates, as do the well known adjusted R2, Mallows Cp, and the AIC criteria. Unfortunately, these corrections do not always work properly. From the objective Bayesian viewpoint, which is mainly based on the results of Section 5.2, two procedures for model selection have been proposed. The main difference between these procedures relies on the form of encompassing – that is, of nesting – the class of all possible submodels. The so called from above procedure (Casella and Moreno 2005, Gir´ on et al. 2004) is based on comparing all the submodels with the full model, and ordering them according to the posterior probabilities of all submodels using the formulae of Section 5.2. We denote these posterior probabilities by Pfa(Mi|Sn) for any submodel Mi, where Sn=(Bn,n,k,k1). The interpretation of these probabilities is that the model with highest posterior probability represents the most plausible reduction in complexity from the full model, the second highest the second most plausible model, and so on. As these posterior probabilities are monote increasing functions of the Bnstatistic for fixed ancillaries, this means that within each class of models Mjthe ordering provided is the same as the one based on the residual sum of squares SSj. Thus, the Bayesian solution from above is the maximal model having the largest posterior probability of its corresponding model. No need for extra adjustment! The so called from below procedure (Gir´ on et al. 2005b and Moreno and Gir´ on 2005), based on the simple fact that the intercept only model is nested into any other possible model as far as it includes the intercept, produces a possibly different ordering of all the submodels of the full model. The ordering is now based on the Bayes factors resulting from comparing the current submodel with the intercept only model. Further, it turns out that this procedure provides a coherent set of model posterior probabilities on the set of all possible submodels denoted by Pfb(Mi|Sn), and these coherent probabilities are monote increasing functions of the R2statistic as now Bn=1−R2,whichin turn, is also a monote decreasing function of the residual sum of squares SSjof the corresponding submodel. This means, as with the from above Bayesian criterion, that the model chosen by the from below criterion is also the maximal model having the largest Bayes factor or, equivalently, the highest model posterior probability Pfb(Mi|Sn). The main conclusion derived from these comparisons is, first, that the two Bayesian criteria always choose a maximal model, i.e. they are compatible with the best subsets El´ ıas Moreno, F. Javier Gir´ on 23 partial ordering and, second, that they are fully automatic in the sense that no tuning of extra parameters, neither the use of outside information nor additional criteria, is needed. Table 3:Comparison of different variable selection criteria for Hald’s data Models From below From above Adjusted Mallows Pfb(Mi|Sn)Pfa(Mi|Sn)R2R2Cp {x1,x2}0.5466 0.7407 0.9787 0.9744 2.6782 {x1,x4}0.1766 0.5364 0.9725 0.9670 5.4958 {x1,x2,x4}0.0889 0.7231 0.9823 0.9764 3.0182 {x1,x2,x3}0.0879 0.7211 0.9823 0.9764 3.0413 {x1,x3,x4}0.0708 0.6809 0.9813 0.9750 3.4968 {x2,x3,x4}0.0165 0.3780 0.9728 0.9638 7.3375 Table 3 compares the results of several model selection criteria for the famous Hald’s data on the composition of cement. The analysis illustrates the Occam’s razor property of the Bayesian criteria. Note also that, for these data, the adjusted R2does not adjust the ordering provided by the original R2for the most plausible models. A large simulation study, see Moreno and Gir´ on (2005b) for the description and extent of the study, has shown that the adjusted R2performs very poorly in almost all situations, a well known fact. Mallow’s Cpand the AIC criteria perform in a very similar way – another well known fact – but they show a poorer behaviour when compared with either Bayesian criteria in most circumstantes. This suggests that the Bayesian criteria account for the difference in dimensionality in some automatic way, hidden in the formulae of their corresponding Bayes factors, in the same manner they automatically obey Occam’s razor principle. Table 4 illustrates these comments for a medium size linear model, k=6, i.e. with five covariates excluding the intercept, sample size n=40 and values of k1ranging from 2to6. Table 4:Comparison of different variable selection criteria for the simulated data. N.oof covariates Criterion 1 2 3 4 5 From below 0.901 0.910 0.962 0.977 1.000 From above 0.573 0.657 0.793 0.927 1.000 Mallows Cp0.500 0.563 0.692 0.850 1.000 Adjusted R20.234 0.292 0.452 0.716 1.000 The model considered for simulation is y=Xα+ε where yis a vector of length 40, Xis a 40 ×6 matrix whose entries were obtained by simulation from a standard normal distribution N(0,1), except the entries in the first 24 On the frequentist and Bayesian approaches to hypothesis testing column which were set equal to 1 to include the intercept, and α α α=(α1,α 2,...,α 6)tis a vector of length 6. The error terms in ε ε ε=(ε1,...,ε n) are i.i.d. εi∼N(0,1). After fixing X, samples of size 5000 were simulated from the model for five different settings of the vector of regression coefficients α α αincluding 1, 2 3, 4 and 5 non zero coefficients. In particular, we set α α α1=(−1,−2,0,0,0,0) α α α2=(−1,−2,2,0,0,0) α α α3=(−1,−2,2,−3,0,0) α α α4=(−1,−2,3/2,−2,2,0) α α α5=(−1,−2,3/2,−2,2,−1). The entries in Table 4 represent the proportion of times that the true model is selected in the first place in the 5000 simulations according to the four criteria and to the number of nonzero regression coefficients in the model. The relation between the p-values, or any equivalent model selection procedure, and the Bayesian model posterior probabilities in the variable selection problem can be summarized as follows. For a fixed sample size n, models with the same number of regressors k1are ordered in the same manner by all criteria: p-values, R2and adjusted R2, Mallows Cp, AIC and the two Bayesian procedures. However, when comparing models with different number of regressors all frequentist and Bayesian criteria generally provide different orderings of the models. But, as we have learned from the simulations, the frequentist behaviour of the Bayesian criteria generally outperforms that of the frequentist ones. Comparisons between the from below and from below Bayesian criteria are discussed in Moreno and Gir´ on (2005b). The conclusion in that paper is that for models with a small or medium number of relevant covariates, as the one illustrated in Table 3, the from below criterion performs better that the from above one, but for models with a large number of influential covariates, the opposite may happen. 6 Discussion Two measures of evidence for hypothesis testing, frequentist and Bayesian, have been considered and compared in this paper for some important testing problems. The case of the standard normal model is not dealt with in the paper as it is a particular case of the normal linear model with no covariates; in this case the sample size is the only ancillary, and the frequentist-Bayesian comparison or calibration just depends on the sample size. We have first recalled that standard normal-theory does not apply to the BehrensFisher problem of testing the equality of the means of two normal populations under heteroscedasticity, as there is no clear-cut p-value and, consequently, frequentist theory has to resort to computing an approximate p-value by adjusting the degrees of freedom of a t-distribution on which the solution of the homoscedastic case is based. El´ ıas Moreno, F. Javier Gir´ on 25 For this problem, we have illustrated the fact that p-values reject the null hypothesis for data for which the Bayesian inference accepts, more markedly as both sample sizes increase. To make p-values more unsatisfactory we have also seen that for small sample sizes p-values would accept the null for data for which the Bayesian rejects. Therefore, we have to admit that the usual interpretation of a p-value as a measure of evidence against the null regardless the sample size, though it entails a notable simplification, may produce wrong answers. It is clear that the enormous success of the p-values in the realm of applications is partially due to their simplicity for scientific communication, but the bad news is that such a simplicity may be misleading. The study of the relation between the two measures of evidence has been extended to normal data in the presence of covariates, that is to the normal linear model. Here two new ancillaries, in addition to the sample size, arise: the number of covariates and the dimension of the null. We have illustrated that the dimension of the null hypothesis is another fundamental ancillary to be taken into account when interpreting p-values. The disregard of this ancillary may produce the rejection of null hypotheses that have high posterior probabilities or the acceptation of nulls that have low posterior probabilities. Finally, we have considered the variable selection problem, a challenging multiple testing problem, because for this problem the sample size and the dimension of the null play a very important role. For selecting variables we necessarily have to jump among models whose null parameter spaces have different dimensions. To overcome this difficulty, the frequentist approach has to adjust some statistic, usually the one based on the ratio of sums of the residuals under the full and null models, in several ways to account for the different number of covariates involved. This accommodation, however, is not very convincing. The objective Bayesian approach seems to deal with this problem in a more appropiate way than the frequentist counterpart because all the ancillaries in the problem are properly taken into account. Furthermore, the objective Bayesian solution works in a fully automatic way in the sense that there is no need for adjusting any tuning parameters. References Berger, J.O. (1985). Statistical Decision Theory and Bayesian Analysis. New York: Springer-Verlag. Berger, J.O. (1994). An overview of robust Bayesian analysis (with discussion). Test, 3, 5-124. Berger, J.O. and Bernardo, J.M. (1992). On the development of the reference prior method. In Bayesian Statistics 4. J.M. Bernardo, J.O. Berger, A.P. Dawid and A.F:M. Smith, eds. Oxford University Press: Oxford, pp. 35-60. Berger, J.O. and Mortera, J. (1999). Default Bayes factors for nonnested hypothesis testing. Journal of the American Statistical Association, 94, 542-554. Berger, J.O. and Pericchi, L.R. (1996). The intrinsic Bayes factor for model selection and prediction. Journal of the American Statistical Association, 91, 109-122. Berger, J.O. and Sellke, T. (1987). Testing a point null hypothesis: the irreconcilability of p-values and evidence (with discussion). Journal of the American Statistical Association, 82, 112-139. 26 On the frequentist and Bayesian approaches to hypothesis testing Bernardo, J.M. and Smith, A.F.M. (1994). Bayesian Theory, New York: Wiley. Box, G.E.P. and Tiao, G.C. (1973). Bayesian Inference in Statistical Analysis. Reading, MA: AddisonWesley. Cano, J.A., Kessler, M. and Moreno, E. (2004). On intrinsic priors for nonnested models. Test, 13, 445463. Casella, G. and Berger, R.L. (1987). Reconciling Bayesian and frequentist evidence in the one-sided testing problem. Journal of the American Statistical Association, 82, 106-111. Casella, G. and Berger, R.L. (1990). Statistical Inference. Belmont: Wadsworth. Casella, G. and Moreno, E. (2005). Intrinsic meta analysis of contingency tables, Statistic in Medicine, 24, 583-604. Casella, G. and Moreno, E. (2005). Objective Bayesian variable selection. Journal of the American Statistical Association (to appear). Fisher, R.A. (1936). The fiducial arguments in statistical inference. Annals of Eugenics, 6, 391-422. Gir´ on, F.J., Mart´ ınez, L. and Imlahi, L. (1999). A characterization of the Behrens-Fisher distribution with applications to Bayesian inference. Comptes rendus de l’Acad´emie des sciences de Paris, Ser I, 701706. Gir´ on, F., Mart´ ınez, L., Moreno, E. and Torres, F. (2003). Bayesian analysis of matched pairs in the presence of covariates. In Bayesian Statistics 7.J.M.Bernardoet al. (eds.), 553-563, Oxford: Oxford University Press. Gir´ on, F.J., Mart´ ınez, L., and Moreno, E. (2003). Bayesian analysis of matched pairs. Journal of Statistical Planning and Inference, 113, 49-66. Gir´ on, F.J., Mart´ ınez, M.L., Moreno, E. and Torres, F. (2004). Objective testing procedures in linear models. Calibration of the p-values. Submitted. Gir´ on, F.J., Moreno, E. and Mart´ ınez, L. (2005). An objective Bayesian procedure for variable selection in regression. In Advances on distribution theory, order statistics and inference.Eds.N.Balakrishnan et al., Birkhauser Boston, (to appear). Jeffreys, H. (1961). Theory of Probability. London: Oxford University Press. Lempers, F.B. (1971). Posterior Probabilities of Alternative Linear Models. Rotterdam: Rotterdam University Press. Lehman, E. (1986). Testing Statistical Hypotheses, (second edition). New York: Wiley. Lindley, D.V. (1970). An Introduction to Probability and Statistics from a Bayesian Viewpoint (Vol. 2). Cambridge: Cambridge University Press. Miller, A.J. (2002). Subset Selection in Regression. 2nd edition. London: Chapman and Hall. Moreno, E. (1997). Bayes Factor for Intrinsic and Fractional Priors in Nested Models: Bayesian Robustness. IMS Lectures Notes-Monograph Series, 31, 257-270. Moreno, E. (2005). Objective Bayesian analysis for one-sided testing, Test, 14, 181-198. Moreno, E., Bertolino, F. and Racugno, W. (1998). An intrinsic limiting procedure for model selection and hypotheses testing. Journal of the American Statistical Association, 93, 1451-1460. Moreno, E., Bertolino, F. and Racugno, W. (1999). Default Bayesian analysis of the Behrens-Fisher problem. Journal of Statistical Planning and Inference, 81, 323-333. Moreno, E., Bertolino, F. and Racugno, W. (2000). Bayesian model selection approach to analysis of variance under heterocedasticity, Journal of the Royal Statistical Society, Series D (The Statistician), 49, 1-15. Moreno, E., Bertolino, F. and Racugno, W. (2003). Bayesian inference under partial prior information, Scandinavian Journal of Statistics, 30, 565-580. Moreno, E. and Cano J.A. (1989). Testing a point null hypothesis: asymptotic robust Bayesian analysis with respect to the priors given on a subsigma field. International Statistical Review, 57, 221-232. El´ ıas Moreno, F. Javier Gir´ on 27 Moreno, E. and Gir´ on, F.J. (2005a). Consistency of Bayes factors for linear models, Comptes rendus de l’Acad´emie des sciences de Paris, Ser I, 911-914. Moreno, E. and Gir´ on, F.J. (2005b). Comparison of Bayesian objective procedures for variable selection in regression. Submitted. Moreno, E., Gir´ on, F.J., and Torres, F. (2003). Intrinsic priors for hypothesis testing in normal regression models. Revista de la Real Academia de Ciencias Serie A,. Mat., 97, 53-61. Moreno, E. and Liseo, B. (2003). A default Bayesian test for the number of components of a mixture. Journal of Statistical Planning and Inference, 111, 129-142. Moreno, E., Torres, F., and Casella, G. (2005). Testing the equality of regression coefficients in heteroscedastic normal regression models. Journal of Statistical Planning and Inference, 131, 117134. O’Hagan, A. (1995). Fractional Bayes factor for model comparison (with discussion). Journal of the Royal Statistical Society Series B, 57, 99-138. O’Hagan, A. and Forster, J. (2004). Bayesian Inference. Kendall’s Advances Theory of Statistics (Vol. 2B). London: Arnold. Pawitan, Y. (2001). In All Likelihood: Statistical Modelling and Inference Using Likelihood. Oxford University Press: Oxford. San Martini, A. and Spezaferri, F. (1984). A predictive model selection criterion. Journal of the Royal Statistical Sociery, Series B, 46, 296-303. Sellke, T., Bayarri, M.J. and Berger, J. (2001). Calibration of p-values for testing precise null hypotheses. The American Statistician, 55, 62-71. Wald, A. (1955). Testing the difference between the means of two normal populations with unknown standard deviations. In Selected papers in Statistics and Probability, T.W. Anderson et al.(eds.), 669-695, Stanford University Press. Wilks, S.S. (1962). Mathematical Statistics. New York: Wiley Welch, B.L. (1951). On the comparison of several means values: an alternative approach. Biometrika, 38, 330-336. Zellner, A. (1986). On assessing prior distributions and Bayesian regression analysis with g-prior distributions. In Bayesian Inference and Decision Techniques: Essays in Honour of Bruno de Finetti, P.K. Goel and A. Zellner (eds.), 233-243, Amsterdam: Elsevier.