Goodness of fit for the logistic regression model using relative belief.

A logistic regression model is a specialized model for product-binomial data. When a proper, noninformative prior is placed on the unrestricted model for the product-binomial model, the hypothesis H 0 of a logistic regression model holding can then be assessed by comparing the concentration of the posterior distribution about H 0 with the concentration of the prior about H 0. This comparison is effected via a relative belief ratio, a measure of the evidence that H 0 is true, together with a measure of the strength of the evidence that H 0 is either true or false. This gives an effective goodness of fit test for logistic regression.

Introduction

Suppose there is a response Y∈{0,1} related to k predictors (X 1,…,X ) via the logistic regression model p(x ′ β)=P(Y=1 | X 1=x 1,…,X =x ) where

and x=(x 1,…,x ), β=(β 1,…,β )∈R . While the use of this model is quite common, the question concerning whether or not the model actually holds has not been fully dealt with in the literature. It is our purpose here to develop a Bayesian approach to this problem.

It is to be noted that, irrespective of whether or not (1) holds, Y | X 1=x 1,…,X =x ∼ Bernoulli (θ(x)) for some θ(x)∈ [0,1] and, if a sample of n(x) is taken at these settings of the predictors, then binomial (n(x),θ(x)). That such data is indeed binomial (independence and constant probability) can be assessed via a runs test applied to each such subsample when n(x)>1. With random sampling from a large enough population the binomial assumption is surely approximately correct and so this aspect of possible model failure is ignored here. The question of interest is whether or not θ(x) is given by (1), at least to a reasonable approximation, and this is the logistic regression assumption.

When there are only categorical predictors, so the X correspond to dummy variables, then indeed (1) holds as then only one of the x =1 with the rest equal to 0 and the relevant probability is exp{β }/(1+ exp{β }) where β is the log of the odds in favor of 1. So in this case the logistic regression is just a reparameterization of the product-binomial model. Goodness of fit beyond the runs tests is then not relevant.

The case when there is at least some quantitative predictors is thus the one of interest. It is to be noted that in any well-designed study there should always be replication, namely, n(x)>1 for some of the x , for precisely the reason that model checking is a necessary part of any statistical analysis in a scientific context. It may be, however, that the data was collected in a somewhat haphazard way, and while model checking is still a requirement, it can’t be expected that this will be as effective as in the designed context.

The general approach taken here can be described as follows. First a noninformative prior is placed on the θ(x). A specific meaning is applied to the word noninformative here, namely, we require that there is no possibility of there being prior-data conflict. Prior-data conflict loosely means that the prior places the bulk of its mass in a region of the parameter space which the data identifies as being unlikely to contain the true value. Prior-data conflict can be measured in a number of ways and this is discussed in Section 2.

Suppose that there are m values for x so the full space for the θ(x) is [0,1]. Let H 0⊂ [0,1] be the subset that corresponds to θ(x)=p(x ′ β) for some β∈R . The prior on the θ(x) leads to a posterior distribution for these quantities. Intuitively, if the posterior is more concentrated about H 0 than the prior, then this is evidence in favor of H 0 with the opposite holding when the posterior is less concentrated about H 0 than the prior. Once a method of measuring concentration about H 0 is selected, this evidence can be measured via a relative belief ratio. While a relative belief ratio is somewhat like a Bayes factor, it will be seen to differ in some key ways. Furthermore, a measure of the strength of this evidence, whether for or against H 0, is presented. Several natural measures of concentration are considered. This is all discussed in Section 3. In Section 4 aspects of the computations are considered including application to a number of examples.

Tsutukawa and Lin (1986) and Bedrick et al. (1996, 1997) are concerned with the Bayesian analysis of logistic regression models although not with goodness of fit. With a p-value based on asympotics, a commonly used goodness of fit statistic for logistic regression is the deviance statistic which is twice the difference between the maximized log-likelihood with no constraints and the maximized log-likelihood, assuming the logistic regression model holds. Chen and Chen (2004) also propose a frequentist asymptotic goodness of fit test in the context of case-control studies. It is shown here that a Bayesian goodness of fit test arises very naturally and that it has a number of advantages. In particular, evidence can be obtained in favor of the logistic regression model, as opposed to only evidence against as with p-values, and there is no appeal to asymptotics.

The Model and prior

Supposed there are m observations s(X)=(s(x 1),…,s(x )) where s(x )∼ binomial (n(x ),θ(x )), the s(x ) are independent, θ(X)=(θ(x 1),…,θ(x ))′∈[0,1] and X=(x 1,…,x )′∈R . Intuitively, a noninformative prior for the θ(X ) is then given by the uniform distribution on [0,1] as this allows for any of the possible values and they are all weighted equally. Of course, other definitions can be provided for noninformativity but this definition seems most suitable for this context as it possesses a key property.

To see this suppose that we have a statistical model {f :θ∈Θ} for the data x, and a prior π. Let T denote a minimal sufficient statistic for the model with marginal model {f :θ∈Θ}. Then is the prior predictive density of T. In Evans and Moshonov (2006) a basic check on the prior is to compute the tail probability M (m (t)≤m (T(x))) and conclude that prior-data conflict exists whenever this probability is small as this implies that the observed data is a priori surprising. The consistency of this check, under quite general conditions, is established in Evans and Jang (2011). If, for example, M (m (t)≤m (T(x)))≡1 for all possible T(x), then it is never the case that there is prior-data conflict and the prior can be called noninformative. Clearly this criterion for noninformativity can be weakened but this is all that is required here.

For the product-binomial, T can be taken equal to s(X). Since the counts are independent and the priors on the θ(x ) are independent and uniform, the prior predictives of the s(x ) are independent. The prior predictive for s(x ) is easily seen to equal 1/(n(x )+1), namely, it is uniform on {0,1,…,n(x )}. As such, for the product-binomial with a uniform prior, M (m (t)≤m (T(x)))≡1 and our criterion for noninformativity is satisfied. The posterior on θ(X) induced by this prior and the observed data s(X) is a product of beta distributions where θ(x ) | s(x )∼ beta (s(x )+1,n(x )−s(x )+1). Note that is easy to generate from both the prior and posterior of θ(X).

Hypothesis assessment via relative belief and concentration

Suppose for the model {f :θ∈Θ} it is desired to assess the hypothesis H 0⊂Θ. In many cases there is a parameter of interest ψ=Ψ(θ) with H 0=Ψ −1{ψ 0}. The evidence for or against H 0 can then be assessed via a relative belief ratio R B (ψ 0 | x)=π (ψ 0 | x)/π (ψ 0), the ratio of the posterior to prior densities of Ψ evaluated at the hypothesized value. If R B (ψ 0 | x)>1, then there is evidence in favor of H 0, as the posterior probability for ψ 0 is greater than the prior probability for ψ 0, while R B (ψ 0 | x)<1 implies there is evidence against H 0 and when R B (ψ 0 | x)=1 there isn’t evidence either way. The strength of the evidence is measured by the posterior probability

For when R B (ψ 0 | x)<1 and (2) is small, then there is a strong belief that the true value of Ψ has a larger relative belief ratio than ψ 0 and so the evidence against ψ 0 is strong. When R B (ψ 0 | x)<1 and (2) is large, then there is only weak evidence against H 0 as this says that there is a large belief that the true value of ψ has the value of its relative belief ratio no greater than R B (ψ 0 | x). When R B (ψ 0 | x)>1 and (2) is large, then there is a weak belief that the true value of Ψ has a larger relative belief ratio than ψ 0 and so the evidence in favor of ψ 0 is strong. Note that in the set {ψ:R B (ψ | x)≤R B (ψ 0 | x)} the value ψ 0 has the most evidence in its favor when R B (ψ 0 | x)>1. When R B (ψ 0 | x)>1 and (2) is small, then there is only weak evidence in favor H 0 as this says that there is a large belief that the true value of ψ has the value of its relative belief ratio greater than R B (ψ 0 | x). The relative belief ratio is discussed in Baskurt and Evans (2013) and a full development of a theory of inference based on this is presented in Evans (2015).

The interpretation of the relative belief ratio as the evidence demands that a relative belief ratio greater than 1 be interpreted as evidence in favour of the hypothesis no matter how much greater it is than 1, assuming it is computed exactly. This is because R B (ψ 0 | x)>1 occurs only when the posterior probability is greater than the prior probability of the hypothesis and that is the basic criterion for saying the data has led to evidence in favor. A similar comment applies to evidence against, namely, when R B (ψ 0 | x)<1. To understand what R B (ψ 0 | x)=1 means consider a discrete context as then this occurs iff the posterior probability of {ψ 0} equals the prior probability of {ψ 0} and this occurs iff the events {ψ 0} and {x} are statistically independent in the joint probability model for (θ,x). In other words R B (ψ 0 | x)=1 iff the actual observed data x tells us nothing about the hypothesis that the true value of ψ=Ψ(θ) is ψ 0. This is clearly a very unusual circumstance but one can easily construct such situations generally when the model contains nonidentifiability. In essence the relative belief ratio is giving the correct assessment of evidence in such a case.

The size of R B (ψ 0 | x) does not necessarily reflect the strength of the evidence. Note that in the discrete case R B (ψ 0 | x)=π (ψ 0 | x)/π (ψ 0)≤1/π (ψ 0) so there is an upper bound on this value. From this it is seen that relative belief ratios do not measure evidence on an absolute scale but they need to be calibrated in each context and that is the role of the strength (2). So even if R B (ψ 0 | x)=1.000005 the strength could be high when (2) is close to 1 as this says our belief that the true value of ψ has a larger relative belief ratio is small. Note that the strength is playing the role of the standard error here as it measures how reliable we believe our assessment of the evidence is. Also, it can happen that even though R B (ψ 0 | x) is very high, (2) can be very small and so the evidence is only weak evidence in favor. This phenomenon is associated with the Jeffreys-Lindley paradox as is discussed in Evans (2015). In short, relative belief ratios need to be calibrated and the calibration depends on the context. The issues concerning measuring strength are a somewhat more involved than the measure of the evidence itself and additional discussion can be found in Evans (2015).

In a number of situations H 0 does not arise via H 0=Ψ −1{ψ 0} for some Ψ in an obvious way and also Π(H 0)=0. The prior nullity may arise because H 0 is a lower dimensional subset of Θ and not because there is no belief that H 0 is true. This is the case with logistic regression when k1]. In such a context it is reasonable to choose where is a measure of the distance of θ(X) to H 0. So with ψ 0=0 then H 0=Ψ −1{ψ 0} and H 0 can be assessed using relative belief. Note that it is clear that in assessing H 0 a comparison is being made between the concentrations of the prior and posterior about H 0. If R B (ψ 0 | x)>1, then the data has lead to the posterior being more concentrated about H 0 than the prior. If R B (ψ 0 | x)<1, then the data has lead to the posterior being less concentrated about H 0 than the prior. The method of concentration, with equal to squared Euclidean distance as discussed in Example 1, was developed for some specific inference problems in Evans et al. (1993, 1997).

While there are many possible choices for two are considered here.

Example 1

Squared Euclidean distance. Let μ(x)= log(θ(x)/(1−θ(x))) denote the logit associated with x and note that the logistic regression model holds iff μ(x)=x ′ β for some β∈R for every x∈R . The logistic regression model thus implies that μ(X)=X β for some β∈R . If a probability distribution is placed on θ(X), then this induces a probability distribution on μ(X) which in turn induces a probability distribution on where it is assumed that X is of full rank. The reason for dividing by the dimension m of μ(X) will become apparent in Example 3 although this clearly has no effect on the optimization. Note that iff the logistic regression model holds for the observed s(X ). So it is natural to measure the concentration of the probability distribution placed on θ(X) about H 0 by seeing how concentrated the induced distribution on is about 0.

Example 2

Kullback-Leibler (KL) distance. The KL distance between the Bernoulli (θ)and the Bernoulli (p)distribution is given by K L(θ,p)=θ log(θ/p)+(1−θ) log((1−θ)/(1−p)). It follows that the KL distance between the product-Bernoulli (θ 1,…,θ ) and the product-Bernoulli (p 1,…,p ) equals So we put and the reason for dividing by m is discussed in Example 4. Given an efficient algorithm for finding the optimal β , it is straightforward to generate from the prior and posterior distributions of

After a choice is made of the methodology proceeds via simulation from the prior and posterior distributions of , then computing and its strength. A significant aspect of this computation is that the prior and posterior densities of typically both vanish at 0 and so the ratio cannot be directly computed. This is a support measure issue arising due to continuity and is dealt with theoretically by defining the relative belief ratio at a point as the limiting ratio of the posterior to prior probabilities of shrinking neighborhoods. Practically this is dealt with by the choice of δ since implies that H 0 holds to the accuracy required in the application. In other words, H 0 holds whenever the difference between the true model and the logistic regression model is of no practical consequence as measured by . The range of the prior distribution of is then discretized via the partition {[0,δ),[δ,2δ),…,[(k−1)δ,k δ)} where k is chosen so that the effective range of this prior distribution is covered. The prior and posterior probability contents of these intervals are estimated by generating large samples from the prior and the posterior distributions of θ, computing for each sampled value, which gives samples from the prior and posterior distributions of and then using the approximate contents for the approximation of the relative belief ratios of the intervals. The relevant relative belief ratio for assessing H 0 is then and the strength of this evidence is assessed by comparing this relative belief ratio against the other values by computing the posterior probability that where the Monte Carlo estimates were used for these computations. It is proved in Evans (2015) that this procedure is consistent as the amount of data increases in the sense that the relative belief ratio converges to the maximum possible value (always greater than 1) and the strength converges to 1 when H 0 is true, and the relative belief ratio converges to 0 and the strength converges to 0 when H 0 is false.

The choice of δ is application dependent as it represents the deviation from the precise null that is just of practical consequence. The two distance measures considered here lead to very natural choices for δ.

Example 3

Squared Euclidean distance and absolute error. For this distance measure let δ equal the maximum squared distance between two logits such that any difference smaller than δ is practically speaking immaterial. In other words, if then this difference is irrelevant from the point of view of the application. It is clear then that for i=1,…,m implies ||μ(X)−X β||2/m < δ while ||μ(X)−X β||2/m < δ implies that the average squared absolute error between individual logits is less than δ. So in practice we proceed by selecting δ and discretizing the prior and posterior distributions of as previously described. Note that it is also reasonable to choose a discretization parameter δ ∗<δ for For example, if δ ∗=δ/m then ||μ(X)−X β||2/m< δ ∗ implies but this might be deemed overly rigorous.

While the interpretation of error in the value of x ′ β is straightforward in linear regression, this is more difficult in logistic regression and it then seems clearer to state bounds on the probabilities. Note, however, for probabilities θ and p, the logits satisfy (log(θ/(1−θ))− log(p/(1−p)))2=(log(θ(1−p)/p(1−θ)))2<δ iff exp(−δ 1/2)<1+(θ−p)/p(1−θ)< exp(δ 1/2), and using e ≈1+x for small x, this is approximately equivalent to (θ−p)2