A default prior distribution for contingency tables with dependent factor levels.

A default prior distribution is proposed for the Bayesian analysis of contingency tables. The prior is specified to allow for dependence between levels of the factors. Different dependence structures are considered, including conditional autoregressive and distance correlation structures. To demonstrate the prior distribution, a dataset is considered which involves estimating the number of injecting drug users in the eleven National Health Service board regions of Scotland using an incomplete contingency table where the dependence structure relates to geographical regions.

Introduction

Contingency tables (e.g.  [1]) are formed when a population is cross-classified according to a series of categories (or factors). Each cell count of the table gives the number observed under each cross-classification. The aim of forming such a table is to summarise the data, and typically, with a view to identifying interactions or relationships between the factors.

The standard statistical practice to model such interactions is the log–linear model (e.g.  [1, Chapter 7]). In this case the logarithm of the expected cell count is proportional to a linear predictor depending on the main effect terms and interaction terms between the factors. Each combination of interaction terms defines its own log–linear model so that the identification of the non-zero interaction terms translates to an exercise in model comparison. Additionally incomplete contingency tables with missing cell counts can be used to estimate closed populations  [4] where some of the factors correspond to sources that have either observed or not observed individuals in the population.

In this paper, we consider the case where the levels of one or more of the factors may be dependent on one another. An obvious example is when one of the factors has levels corresponding to geographical regions or locations which may be dependent due to their geographical proximity. In these cases, we may expect the parameters of the log–linear model to have some dependence structure. Bayesian analysis of contingency tables is common (e.g.  [3], [13], [5]) and is the approach taken here. One feature of the Bayesian approach is that prior information on the interaction terms can be incorporated through the prior distribution. We take the position of having weak prior information on the magnitude of the log–linear parameters but wish to incorporate the information provided by the dependence structure mentioned above. In the case of weak prior information and model uncertainty, care must be taken when specifying prior distributions due to Lindley’s paradox (e.g.  [16, pp. 77–79]). There have been several attempts in the literature (e.g.  [3], [15], [18]) to specify “default” prior distributions that can be applied for log–linear models under model uncertainty. We extend these approaches by developing a default prior that can take account of the dependence structure between the factor levels and can be seen as a generalisation of the above mentioned priors. The proposed prior is constructed by conditioning on the constraints on the parameters which are introduced in contingency table analysis to maintain identifiability of the parameters.

This paper is organised as follows. In Section  2 we set out our notation and briefly describe log–linear models. In Section  3 we derive our proposed default prior distribution including descriptions of different dependence structures. Finally, we apply our proposed prior to a real data application in Section  4, which involves estimating the number of injecting drug users in Scotland. Here, one of the factors corresponds to geographical regions, and we wish to take account of the possible dependence structure that may exist for the regions.

Notation and log–linear models

Notation

We assume that there are a total of factors such that each factor has levels. The corresponding contingency table has cells. Let be the vector of cell counts with elements denoted as and where identifies the combination of factor levels that cross-classify the cell . Let be set of all cross-classifications so that Finally, let be the total population size. In the case of an incomplete contingency table, is unknown, since elements of are unknown.

As a pedagogic example that we use for illustrative purposes throughout, suppose that there are three factors used to cross-classify a population of hospital patients: age (2 levels: young; old), hypertension (2 levels: no; yes) and region (3 levels: A; B; C). In this example, , where and , and the three factors (age, hypertension and region) have been labelled 1, 2 and 3, respectively. It follows that there are cells.

Log–linear models

We now briefly describe log–linear models and initially assume that the form of the log–linear model is known, i.e. it is known which interactions are present. We extend to the case of model uncertainty later in this section. Let denote the linear predictor associated with cell , where with denoting the intercept term, the vector of log–linear parameters (i.e. the main effects and interaction terms) and the vector of zeros and ones identifying which elements of are applicable to cell .

For identifiability, certain elements of are constrained, e.g. by sum-to-zero, or corner-point constraints, so we can rewrite as where is the vector of unconstrained regression parameters, and is the vector which identifies which elements of correspond to cell , with .

Finally, let be the vector with elements , and let be the model matrix with rows . Then we can write where denotes the vector of ones.

For the statistical analysis of contingency tables, it is common to assume that independently, where .

In practice, we typically do not know the form of the log–linear model. This is equivalent to not knowing the elements of and , or the columns of . Let be the set of competing log–linear models which are indexed by . Associated with each log–linear model are and , where and are vectors, and are vectors, and is an matrix.

In the next section, we derive a default prior distribution for . For the intercept, , we assume a prior given by . Although this prior is improper, the resulting posterior is still proper  [5]. This prior will not cause a problem under Lindley’s paradox since it is present for all models in [16, p. 174].

A default prior distribution for

Derivation

In this section we develop a default prior distribution for . For notational simplicity, we drop the dependency on the model by removing the superscript .

Suppose that there are a total of log–linear terms and where , for , is the vector corresponding to the regression parameters for the main effect or interaction term . Similarly let denote the corresponding vector of log–linear parameters, for .

Let be the set of main effect terms that define the -way interaction . Dellaportas and Forster  [3] refer to as the constituent terms of the interaction. Note that and if corresponds to a main effect then has only one element, i.e. . Consider the pedagogic example, from Section  2.1, and corresponding to the 2-way interaction between age and region so that and . The constituent terms, , have two elements: the terms corresponding to age and region.

We initially consider deriving the default prior distribution under sum-to-zero constraints. We describe how the prior can be extended to any system of constraints in Section  3.4. Following Dellaportas and Forster  [3] we assume that has a multivariate normal distribution with mean zero, where and are independent for and . Thus, all that remains is to specify the covariance matrix for each , for .

The elements of are subject to constraints and can be written in the form where is a matrix defining the constraints. Under sum-to-zero constraints, can be written as where is the identity matrix, is a matrix and is a permutation matrix. For corresponding to the age and region interaction in the pedagogic example, The elements of are ordered so that the factor levels of region vary the fastest.

Initially, ignoring the constraints that are applied to , we assume that the distribution of is where and is a positive-definite scale matrix. The off-diagonal elements of control the dependence structure or correlation between the elements of the constrained parameters, , corresponding to different factor levels.

It follows from (2), (3) that Let be the permuted elements of according to the inverse permutation , so that and . The prior distribution for is the conditional distribution of (which is ) given that , i.e. we find the distribution of from (4) subject to the constraints. It can be shown (see Appendix A) that where In the next two sections we consider . It may be that is completely specified a priori. The most plausible situation for this is when we assume independence between the levels of this term and . We consider this case in Section  3.2. In Section  3.3 we also consider where is unknown due to its dependence on some unknown hyperparameter which controls the strength of correlation between the elements of .

Independent correlation structure

Suppose we assume that the factor levels are independent, i.e. , so that Denote by the matrix formed by the columns of corresponding to . Since is a permutation of the matrix formed by stacking to form an matrix, it follows that and therefore . The corresponding prior distribution for is If , then since (under sum-to-zero constraints) , for all   [14], it follows that the prior distribution for is If is unknown and given a prior distribution, then (7) is a hierarchical prior distribution that is identical to the generalised hyper-g prior proposed by Sabanes-Bové and Held  [18] for generalised linear models (GLMs) when applied to log–linear models. If, instead, is fixed then (7) is the default prior distribution considered by Dellaportas and Forster  [3] who advocate setting for some constant , which represents the number of units of prior information. Ntzoufras et al.  [15] use under their unit information prior for GLMs when applied to log–linear models.

General correlation structure

We now consider terms, , whose constituent terms, , contain factors with correlated levels and depends on some unknown hyperparameter . This hyperparameter, , controls the strength of correlation through some structure imposed on . Initially consider a main effect term . In this paper we focus on the case where the factor levels correspond to geographical regions or locations and propose two structural forms for . However there exist many possible applications with correlated factor levels and other correlation structures that can be used depending on the nature of the factor levels.

Conditional autoregressive structure

Suppose that the levels correspond to regions. Let be the neighbourhood matrix with th element for . Then for the conditional autoregressive (CAR) structure (e.g.  [2]), where determines the strength of spatial correlation for the constrained parameters. To ensure that is positive-definite, the hyperparameter must lie in the interval , where and are the maximum and minimum eigenvalues of , respectively.

Distance correlation structure

Suppose the levels correspond to locations such as cities. Then the th element of is given by a correlation function that depends on the distance, , between locations and , and . For example, the Gaussian correlation function gives where, again, controls the strength of correlation.

Note that in both examples, the hyperparameter, , is not actually a correlation coefficient; it merely controls the strength of correlation. We need to specify a prior distribution for . This will depend on the application.

For a term that corresponds to an interaction term, we propose The form given by (8) has been chosen for its consistency. Suppose that the correlation between two levels of a main effect term is . Then, for an interaction involving this main effect, the correlation between the two levels will be if and only if the factor levels of the other constituent terms are identical. To demonstrate this we return to our pedagogic example where the regions A and B, and B and C are neighbours, but A and C are not neighbours. A CAR structure is specified. In this example, the neighbourhood matrix is so that for the main effect of region is The eigenvalues of are , so, therefore, . If an independent correlation structure is specified for the main effect of age, then The correlation between A and B for the main effect of region is . For the age and region interaction, the correlation between levels involving A and B is if and only if they have the same level for age. It now follows from (6), (9) that the scale matrix for the prior distribution is If we denote the regression parameters for this term as , where , then the prior correlation between and is If , corresponding to independence between the regions, i.e. , and thus we have the Sabanes-Bové and Held  [18] prior, then . The function is increasing in but the correlation is always negative. This is caused by the sum-to-zero constraints. As increases, the magnitude of the negative correlation decreases.

A further advantage of using the structure defined by (8) is computational. If we assume that the independence model, containing only the main effect terms, is the simplest model we wish to consider then we will always have the same set of hyperparameters in each model.

Alternative constraint systems

We now consider alternative constraint systems to sum-to-zero constraints, e.g. corner-point or Helmert constraints. Let and denote the vectors of regression parameters under the alternative and sum-to-zero constraints, respectively. Since, under the sum-to-zero constraints, each component, , of has a normal distribution, then has a normal distribution with mean zero and variance matrix . It can be shown (see Appendix B) that where and are the model matrices under the alternative and sum-to-zero constraints, respectively, is the matrix of ones and Therefore , where the prior variance matrix, , is given by Note that, under the alternative constraints, and may no longer, necessarily, be independent. This is equivalent to the fact that (given by the above expression) may no longer, necessarily, be block diagonal.

Under the independence structure described in Section  3.2, where , for , then The matrix is called the hat matrix and is invariant to the type of constraint system used, i.e. and therefore Therefore the proposed prior distribution is a generalisation of the default prior distribution of Sabanes-Bové and Held  [18] for any type of constraint system.

Example: estimating the number of injecting drug users (IDUs) in Scotland from capture–recapture data

In this section we apply our proposed default prior distribution to an incomplete contingency table which has six factors and 352 cells that involves estimating the number of injecting drug users (IDUs) in Scotland in 2006. These data have been previously analysed by King et al.  [12] and Overstall et al.  [17]. The six factors are social enquiry reports (2 levels: observed; unobserved); hospital records (2 levels: observed; unobserved); Scottish drug misuse database (2 levels: observed; unobserved); age (2 levels: ≤35 years; >35 years); gender (2 levels: male; female) and region (11 levels: National Health Service (NHS) board regions—see Fig. 1). The first three factors are sources and the 44 cells which correspond to not being observed by any of these sources for the different age/gender/region combinations have missing counts. Therefore the total population of IDUs, , is unknown. We use Markov chain Monte Carlo (MCMC) methods to obtain posterior distributions for the missing cell entries and therefore a posterior distribution for the total population of IDUs.

Fig. 1

Map showing the eleven regions of Scotland which correspond to National Health Service (NHS) board regions.

King et al.  [12] and Overstall et al.  [17] merged the eleven regions into just two levels: Greater Glasgow and Clyde, and the Rest of Scotland. Without merging, using all eleven distinct regions, there are small cell counts for many of the regions. For instance, in one region there are only 19 observed IDUs over all source, age and gender cross-classifications. This suggests that a prior distribution that involves smoothing (or borrowing of information), such as the prior proposed in Section  3, is required. We apply the proposed prior where the independence structure is specified for all of the factors except region where we use the CAR structure described in Section  3.3. By calculating the eigenvalues of the neighbourhood matrix, , for this example, and . We place a uniform prior on in the interval . The prior distribution for each is where is given by (6). Following from Section  3.2, we set , with where IG denotes the inverse-gamma distribution, and , as suggested by Sabanes-Bové and Held  [18]. We only specify non-zero prior model probabilities for the log–linear models that contain at most two-way interactions and assume a discrete uniform prior over all of these models. It was found that this allowed enough complexity to obtain an adequate overall model when using the Bayesian -value to assess model adequacy (see,  [8, Chapter 6]).

We use the data-augmentation MCMC approach proposed by King and Brooks  [13] with the reversible jump implementation for GLMs of Forster et al.  [6] to make moves between log–linear models and the weighted least squares Metropolis–Hastings implementation of Gamerman  [7] to make moves within the same log–linear model. We ran the algorithm for one million iterations (discarding the first 10% as burn-in).

For the total population size of IDUs, we obtain a posterior distribution for the total population size with a mean of 21 700 and a 95% highest posterior density interval (HPDI) of (18 900, 24 800). Overstall et al.  [17] obtained a posterior mean of 24 000 and a 95% HPDI of (19 500, 29 700) and King et al.  [12] a mean of 25 000 with a 95% HPDI of (20 700, 35 000). The advantage of our approach over the latter two analyses is that we are able to provide posterior distributions of the total population size in each NHS board region, broken down by age and gender. Our approach also results in a smaller credible interval for the total population size due to it allowing for correlated regions and not discarding information by merging the factor levels of region.

The posterior mean of is 0.108 with a 95% HPDI of . The posterior probability of being positive is 0.816. It follows that the Bayes factor in support of the hypothesis that is 8.205. Therefore there appears to be positive evidence  [11] in support of positive spatial correlation between the regions of Scotland.

Concluding remarks

In this paper we have proposed a default prior distribution for the regression parameters of a log–linear model that can take account of any dependence structure that may exist between the factor levels. This prior can be applied in situations of model uncertainty and can be seen as a generalisation of other default prior distributions applied to log–linear models including those of Dellaportas and Forster  [3], Ntzoufras et al.  [15] and Sabanes-Bové and Held  [18].