Generating contingency tables with fixed marginal probabilities and dependence structures described by loglinear models

Abstract

We present a method to generate contingency tables that follow loglinear models with prescribed marginal probabilities and dependence structures. We make use of (loglinear) Poisson regression, where the dependence structures, described using odds ratios, are implemented using an offset term. Other statistical models related to loglinear models that fall into the scope of this paper, such as the logistic regression model, the latent class model and the extended Rasch are discussed as well. We apply this methodology to carry out simulation studies in the context of population size estimation using dual system and triple system estimators, popular in official statistics. These estimators use contingency tables that summarize the counts of elements enumerated or captured within lists that are linked. The simulation is used to investigate these estimators in the situation that the model assumptions are fulfilled, and the situation that the model assumptions are violated.

Keywords:

• Contingency tables
• loglinear model
• odds ratio
• offset
• simulation
• dual-system estimator
• triple-system estimator

1. Introduction

In simulation studies that make use of contingency tables, samples are drawn from a population with properties that are well understood. In this paper a method is proposed for the generation of contingency tables with fixed marginal probabilities and prescribed dependence structures that are defined in terms of odds ratios. An odds ratio describes the strength of the association for a two-way contingency table. For multi-way tables conditional odds ratios play the same role for partial associations between variables [1]. A simple example of such a population would be a two-way contingency table with marginal probabilities of 0.7 and 0.8 and the dependence represented by an odds ratio of 2.

Loglinear models have interaction parameters that can be directly understood in terms of odds ratios. Therefore the proposed method is built around the use of loglinear models. Models with two and three variables are explored, where each variable is dichotomous. However, the methodology presented can generate populations using loglinear models of any size, in terms of both the number of variables and the number of levels of the variables. Also, we can generate populations following logistic regression models having only categorical explanatory variables, and latent variable models that can be formulated as loglinear models, such as the latent class model [2] and the extended Rasch model [3].

There are a few ways to fit loglinear models. It suits our purposes to define loglinear models in terms of Poisson regression models with a log link [4]. Such models can be fitted with, for example, iteratively reweighted least squares (IRLS) [5]. The models we use require the marginal probabilities and odds ratios as input parameters. As a first step, the given marginal probabilities are used to generate the probabilities in a loglinear model for which complete independence holds [6,7, p. 345]. The prespecified interaction structure is imposed by using an offset, which is an additional explanatory variable with a prespecified parameter value that is fixed and not estimated.

When each margin of a contingency table has only two categories, the contingency table is a summary of the observations of correlated binary variates. Previously, Lee [8] proposed a linear programming method for this problem. In comparison to his method, our method is simpler. Also, we present our method in the context of various applications of the loglinear model where conditional odds ratios can also be specified, and of latent variable models. The work of Gange [9] is related to our work in the sense that he uses the iterative proportional fitting algorithm, an algorithm that was very popular in the early ages of loglinear modelling; yet he does not make a link to loglinear models and mentions odds rations only in passing. Choi [10] also presented a method based on loglinear models, where the association is formulated in terms of the Pearson chi-square coefficient, and iterative proportional fitting is used in the generation of probabilities. Thus the way in which the association is specified by Choi is different from our method, that is more closely related to the way in which the loglinear model is defined. Similar to our work, Kateri [11, section 2.2.5, p. 43] described that local odds ratios and row and column marginal distributions uniquely determine the table of joint probabilities, which also holds for any other minimal set of odds ratios. Our work further develops this and demonstrates a simple approach to generating data from this information. Geenens [12] sets out a similar approach into a framework where the dependence is measured using copulas. He explains that the bivariate distributions are uniquely determined by a ‘margin-free’ association parameter (and that multiple different kinds of parameter are possible, but with one-to-one relations between them.

There is also much related but different work on generating correlated binary variates, see, for example, Emrich and Piedmonte [13], Lee [8], Gange [9], Park et al. [14], Lee [15], Lunn and Davies [16], Al Osh and Lee [17], Demirtas [18] and Jiang et al. [19]. There is also some (though rather less) work on generating correlated multinomial or ordinal variables which can be used when the margins have more than two categories, for example Gange [9], Choi [10], Amatya and Demirtas [20] and Touloumis [21]. These approaches mostly focus on generating population probabilities in order to assess the performance of generalized estimating equation (GEE) models. In passing, some of these papers mention that their methodology can also be used with odds ratios instead of correlations, but they do not provide a detailed description of how this works. Geenens [12] extends the connection of copulas and margin-free association measures to general discrete distributions. Combined with marginal information these generate distributions and Geenens demonstrates their existence and uniqueness. Our work is different in that we focus on generating data from loglinear models where the dependence is specified using odds ratios rather than correlations, and describe how this leads to a simple implementation using standard functions for fitting loglinear models.

The methodology we propose is illustrated with an application to population size estimation, as has been used in population censuses. This follows on from the work of Brown et al. [22], Baffour et al. [23] and Gerritse et al. [24] who used a fixed odds ratio to investigate the impacts of dependence on population size estimation. In this application, lists where individuals within the target population are enumerated are used to estimate the population sizes. These lists are subject to incompleteness (undercoverage error), and multiple system estimation has been proposed to adjust for this undercoverage error. Some assumptions about the dependence of the lists must be made (though these may be quite weak where there are many lists), so it is important to test these methods under prespecified conditions. The proposed methodology and its application allow investigations into the impacts of given dependence structures within contingency tables to be carried out, varying the sizes of odds ratios, population sizes and response probabilities, and under different models.

The structure of the paper is as follows. In Section 2 we outline the methodology around log linear models and model fitting. From this we develop the proposed method for including prespecified interaction parameters in Section 2.4. In Section 3 we present an application where the proposed method is used to investigate the properties of population size estimation under dependence. A simulation study is presented for two-way and three-way contingency tables, alongside steps for users to implement the methodology. Section 4 makes an assessment of the advantages of our approach compared with existing methods for correlated binary or categorical data.

2. Methodology

2.1. Loglinear models

2.1.1. Models for two-way tables

The simplest situation that we consider is the loglinear model with two dichotomous variables A and B. Let A be indexed by i = 0, 1 and B by j = 0, 1. We denote the expected count for cell $(i,j)$ of the two-way contingency table of variables A and B by ${\mu }_{\mathrm{ij}}$. We denote the loglinear parameters by λ. Then the saturated loglinear model for two lists is: (1) $\log {\mu }_{\mathrm{ij}}=\lambda +{\lambda }_i^A+{\lambda }_j^B+{\lambda }_{\mathrm{ij}}^{\mathrm{AB}},$(1) where λ is the intercept term, ${\lambda }_i^A$ and ${\lambda }_j^B$ are the respective parameters for variables A and B, and ${\lambda }_{\mathrm{ij}}^{\mathrm{AB}}$ is the two factor interaction parameter. The parameters are identified by setting the parameter equal to 0 when an index is 0, i.e. ${\lambda }_0^A={\lambda }_0^B={\lambda }_{00}^{\mathrm{AB}}={\lambda }_{01}^{\mathrm{AB}}={\lambda }_{10}^{\mathrm{AB}}=0.$ Thus the model has four parameters and four observed counts, and hence is called saturated. The loglinear interaction parameter ${\lambda }_{11}^{\mathrm{AB}}$ is closely related to odds ratio ${\theta }^{\mathrm{AB}}$ and the expected counts by (2) ${\theta }^{\mathrm{AB}}=\text{exp}({\lambda }_{11}^{\mathrm{AB}})=\frac{({\mu }_{11}/{\mu }_{01})}{({\mu }_{10}/{\mu }_{00})}=\frac{{\mu }_{11}{\mu }_{00}}{{\mu }_{01}{\mu }_{10}}.$(2) When ${\lambda }_{\mathrm{ij}}^{\mathrm{AB}}=0$, the odds ratio ${\theta }^{\mathrm{AB}}=1$, and the independence model holds.

The proposed methodology can be extended in a straightforward way to variables with more than two levels (compare [11], section 2). Assume a $3\times 3$ table with levels i, j = 1, 2, 3, with level 3 being the reference category. The model is described in Equation (1(1) $\log {\mu }_{\mathrm{ij}}=\lambda +{\lambda }_i^A+{\lambda }_j^B+{\lambda }_{\mathrm{ij}}^{\mathrm{AB}},$(1) ) with i, j now having different levels than in the $2\times 2$ context. If we identify the 9 interaction parameters by setting these parameters equal to 0 when $(i,j=3)$, then there are four interaction parameters to be estimated, namely ${\lambda }_{11}^{\mathrm{AB}},{\lambda }_{12}^{\mathrm{AB}},{\lambda }_{21}^{\mathrm{AB}}$ and ${\lambda }_{22}^{\mathrm{AB}}$. Thus there are four odds ratios, ${\theta }_{11}^{\mathrm{AB}},{\theta }_{12}^{\mathrm{AB}},{\theta }_{21}^{\mathrm{AB}}$ and ${\theta }_{22}^{\mathrm{AB}}$, where for each cells in level $(i,j=3)$ is involved. For example, the odds ratio (3) ${\theta }_{12}^{\mathrm{AB}}=\exp \left({\lambda }_{12}^{\mathrm{AB}}\right)=\frac{({\mu }_{12}/{\mu }_{13})}{({\mu }_{32}/{\mu }_{33})}=\frac{{\mu }_{12}{\mu }_{33}}{{\mu }_{13}{\mu }_{32}}.$(3)

2.1.2. Models for three-way tables

Consider the loglinear model for three dichotomous variables, where the third variable is denoted as C. Here the saturated model is (4) $\log {\mu }_{\mathrm{ijk}}=\lambda +{\lambda }_i^A+{\lambda }_j^B+{\lambda }_k^C+{\lambda }_{\mathrm{ij}}^{\mathrm{AB}}+{\lambda }_{\mathrm{ik}}^{\mathrm{AC}}+{\lambda }_{\mathrm{jk}}^{\mathrm{BC}}+{\lambda }_{\mathrm{ijk}}^{\mathrm{ABC}},$(4) where λ is the intercept term, ${\lambda }_i^A$ ${\lambda }_j^B$ and ${\lambda }_j^C$ are respectively the parameters for variables A, B and C, ${\lambda }_{\mathrm{ij}}^{\mathrm{AB}}$, ${\lambda }_{\mathrm{ik}}^{\mathrm{AC}}$ and ${\lambda }_{\mathrm{jk}}^{\mathrm{BC}}$ are the two factor interaction parameters and ${\lambda }_{\mathrm{ijk}}^{\mathrm{ABC}}$ is the three factor interaction parameter. In order to identify the parameters, similar identifying restrictions hold as for the two-variable case.

When three dichotomous variables are present, the odds ratio for each pairwise relationship is conditional on the outcome of the variable not included in the pairwise relationship. E.g. when there is a relationship between variables A and B, the odds ratio is conditional on variable C: (5) $\begin{array}{rl}{\theta }^{\mathrm{AB}|C=0}& =\text{exp}({\lambda }_{11}^{\mathrm{AB}})=\frac{{\mu }_{110}{\mu }_{000}}{{\mu }_{100}{\mu }_{010}},\end{array}$(5) (6) $\begin{array}{rl}{\theta }^{\mathrm{AB}|C=1}& =\text{exp}({\lambda }_{11}^{\mathrm{AB}}+{\lambda }_{111}^{\mathrm{ABC}})=\frac{{\mu }_{111}{\mu }_{001}}{{\mu }_{101}{\mu }_{011}}.\end{array}$(6)

2.2. Latent variable models

Loglinear models are also used in the context of latent variable models. The Rasch model, popular in the field of item response theory in psychometrics and in the field of population size estimation, assumes a single continuous latent variable that explains the dependence between observed dichotomous variables, i.e. the observed variables are conditionally independent given the latent variable. A version of the popular Rasch model [25, originally published in 1960] is the so-called extended Rasch model [3], that can be written for three observed variables as model (4(4) $\log {\mu }_{\mathrm{ijk}}=\lambda +{\lambda }_i^A+{\lambda }_j^B+{\lambda }_k^C+{\lambda }_{\mathrm{ij}}^{\mathrm{AB}}+{\lambda }_{\mathrm{ik}}^{\mathrm{AC}}+{\lambda }_{\mathrm{jk}}^{\mathrm{BC}}+{\lambda }_{\mathrm{ijk}}^{\mathrm{ABC}},$(4) ) with all two-factor interactions identical, i.e. ${\lambda }_{11}^{\mathrm{AB}}={\lambda }_{11}^{\mathrm{AC}}={\lambda }_{11}^{\mathrm{BC}}$. In psychometric testing the variables A, B and C stand for item responses, and in the Rasch model individuals are ordered by the difficulty of the items and items are ordered by the aptitude of the individuals. For example, an individual may master item A but neither B nor the most difficult item C. Another individual masters A and B and thus has a higher aptitude than the first individual. Thus the model allows individuals to have heterogeneous probabilities to master items.

The latent class model assumes the existence of a categorical latent variable X with levels x, and given the level x the observed variables are independent. Let's assume that there are four observed variables A, B, C and D (having levels indexed by l) and a dichotomous variable X. Written in terms of (conditional) probabilities we have, (7) ${\pi }_{\mathrm{ijkl}}^{\mathrm{ABCD}}={\mathrm{\Sigma }}_x{\pi }_{\mathrm{ijklx}}^{\mathrm{ABCDX}}={\mathrm{\Sigma }}_x{\pi }_x^X{\pi }_{i|x}^{\mathrm{AX}}{\pi }_{j|x}^{\mathrm{BX}}{\pi }_{j|x}^{\mathrm{CX}}{\pi }_{l|x}^{\mathrm{DX}}.$(7) This model was originally proposed by Lazersfeld and Henry [26]. Haberman [2] showed that the model can also be written as a loglinear model: (8) $\log {\mu }_{\mathrm{ijklx}}^{\mathrm{ABCDX}}=\lambda +{\lambda }_i^A+{\lambda }_j^B+{\lambda }_k^C+{\lambda }_l^D+{\lambda }_{\mathrm{ix}}^{\mathrm{AX}}+{\lambda }_{\mathrm{jx}}^{\mathrm{BX}}+{\lambda }_{\mathrm{kx}}^{\mathrm{CX}}+{\lambda }_{\mathrm{lx}}^{\mathrm{DX}},$(8) with ${\mu }_{\mathrm{ijkl}}^{\mathrm{ABCD}}={\mathrm{\Sigma }}_x{\mu }_{\mathrm{ijklx}}^{\mathrm{ABCDX}}$. For the purpose of this paper, writing the latent class model as in Equation (8(8) $\log {\mu }_{\mathrm{ijklx}}^{\mathrm{ABCDX}}=\lambda +{\lambda }_i^A+{\lambda }_j^B+{\lambda }_k^C+{\lambda }_l^D+{\lambda }_{\mathrm{ix}}^{\mathrm{AX}}+{\lambda }_{\mathrm{jx}}^{\mathrm{BX}}+{\lambda }_{\mathrm{kx}}^{\mathrm{CX}}+{\lambda }_{\mathrm{lx}}^{\mathrm{DX}},$(8) ) may not be that advantageous to generate a population, as Equation (7(7) ${\pi }_{\mathrm{ijkl}}^{\mathrm{ABCD}}={\mathrm{\Sigma }}_x{\pi }_{\mathrm{ijklx}}^{\mathrm{ABCDX}}={\mathrm{\Sigma }}_x{\pi }_x^X{\pi }_{i|x}^{\mathrm{AX}}{\pi }_{j|x}^{\mathrm{BX}}{\pi }_{j|x}^{\mathrm{CX}}{\pi }_{l|x}^{\mathrm{DX}}.$(7) ) can be used for this purpose in a straightforward way by specifying the (conditional) probabilities on the r.h.s. of Equation (7(7) ${\pi }_{\mathrm{ijkl}}^{\mathrm{ABCD}}={\mathrm{\Sigma }}_x{\pi }_{\mathrm{ijklx}}^{\mathrm{ABCDX}}={\mathrm{\Sigma }}_x{\pi }_x^X{\pi }_{i|x}^{\mathrm{AX}}{\pi }_{j|x}^{\mathrm{BX}}{\pi }_{j|x}^{\mathrm{CX}}{\pi }_{l|x}^{\mathrm{DX}}.$(7) ), thus generating population probabilities for the four-way table.

2.3. Logistic regression

Due to the relation between the loglinear model and the logistic regression model with categorical explanatory variables [6, section 2.2.4], the method proposed in this paper can also be used for this subset of logistic regression models. We show the relation between the models for a simple example. Consider model (4(4) $\log {\mu }_{\mathrm{ijk}}=\lambda +{\lambda }_i^A+{\lambda }_j^B+{\lambda }_k^C+{\lambda }_{\mathrm{ij}}^{\mathrm{AB}}+{\lambda }_{\mathrm{ik}}^{\mathrm{AC}}+{\lambda }_{\mathrm{jk}}^{\mathrm{BC}}+{\lambda }_{\mathrm{ijk}}^{\mathrm{ABC}},$(4) ). Let the variable C with levels k = 0, 1 be the response variable and variables A and B be categorical variables. Then the logistic regression model is (9) $\log \frac{{\mu }_{\mathrm{ij}1}}{{\mu }_{\mathrm{ij}0}}=\log {\mu }_{\mathrm{ij}1}-\log {\mu }_{\mathrm{ij}0}={\lambda }_k^C+{\lambda }_{\mathrm{ik}}^{\mathrm{AC}}+{\lambda }_{\mathrm{jk}}^{\mathrm{BC}}+{\lambda }_{\mathrm{ijk}}^{\mathrm{ABC}},$(9) where the λ-parameters ${\lambda }_k^C,{\lambda }_{\mathrm{ik}}^{\mathrm{AC}},{\lambda }_{\mathrm{jk}}^{\mathrm{BC}}$ and ${\lambda }_{\mathrm{ijk}}^{\mathrm{ABC}}$ can be replaced by β-parameters (not involving k and C) $\beta ,{\beta }_i^A,{\beta }_j^B$ and ${\beta }_{\mathrm{ij}}^{\mathrm{AB}}$ respectively. The loglinear model $[\mathrm{AB}][C]$ corresponds with the intercept only model in logistic regression as ${\lambda }_{\mathrm{ik}}^{\mathrm{AC}}={\lambda }_{\mathrm{jk}}^{\mathrm{BC}}={\lambda }_{\mathrm{ijk}}^{\mathrm{ABC}}=0$ implies ${\beta }_i^A={\beta }_j^B={\beta }_{\mathrm{ij}}^{\mathrm{AB}}$ = 0.

Multinomial and ordinal response models can also be dealt with using this approach if they have equivalent loglinear models. Examples include the baseline category-logit model and the continuation-ratio logit model, see Agresti [1].

2.4. The proposed method

Loglinear models are usually estimated using the maximum likelihood method. An important property of maximum likelihood estimates of loglinear models is that, for the parameters fitted, the corresponding margins of the observed counts equal the margins of the fitted counts [6, p. 69]. These fitted counts are the sufficient statistics. We denote maximum likelihood estimates by adding a hat to parameters and expected counts.

We focus here on the estimation of hierarchical loglinear models, starting with the two-way table and extending to three-way tables later on. In hierarchical loglinear models, when higher-order parameters are included, the lower order parameters are included as well, for example, if ${\lambda }_{\mathrm{ij}}^{\mathrm{AB}}$ is included, the lower order parameters ${\lambda }_i^A$ and ${\lambda }_j^B$ are included too. The use of hierarchical loglinear models allows for simple and easier interpretation of the model parameters.

2.4.1. Model fitting for two-way tables

For a two-way contingency table with observed counts $n_{\mathrm{ij}}$, when we fit model (1(1) $\log {\mu }_{\mathrm{ij}}=\lambda +{\lambda }_i^A+{\lambda }_j^B+{\lambda }_{\mathrm{ij}}^{\mathrm{AB}},$(1) ), the fitted values are equal to the observed counts, i.e. ${\hat{\mu }}_{\mathrm{ij}}=n_{\mathrm{ij}}$. This model is denoted by the variables constituting the highest fitted margins, i.e. [AB]. The margins $n_{\mathrm{ij}}$ are the sufficient statistic for the saturated model. For the independence model, denoted by [A][B], the margins of variable A and variable B are fitted, i.e. ${\hat{\mu }}_{\mathrm{ij}}={\hat{\mu }}_{i+}{\hat{\mu }}_{+j}/{\hat{\mu }}_{++}=n_{i+}{n}_{+j}/n_{++}$, where a ‘+’ denotes that the sum is taken over the corresponding index.

Although the two most relevant models, the saturated model and the independence model, can be fitted in a straightforward way, we fit the model here through loglinear Poisson regression, as this suits our purposes later on. Assume that the four elements $n_{\mathrm{ij}}$ are counts derived from four Poisson distributions. Then a likelihood can be set up for ${\mu }_{\mathrm{ij}}$ and maximized over the parameters. If we collect the elements ${\mu }_{\mathrm{ij}}$ into a column vector, then model (1(1) $\log {\mu }_{\mathrm{ij}}=\lambda +{\lambda }_i^A+{\lambda }_j^B+{\lambda }_{\mathrm{ij}}^{\mathrm{AB}},$(1) ) can be written as, (10) $\left[\begin{array}{c}\log ({\mu }_{11})\\ \log ({\mu }_{10})\\ \log ({\mu }_{01})\\ \log ({\mu }_{00})\end{array}\right]=\left[\begin{array}{cccc}1& 1& 1& 1\\ 1& 1& 0& 0\\ 1& 0& 1& 0\\ 1& 0& 0& 0\end{array}\right]\left[\begin{array}{c}\lambda \\ {\lambda }_1^A\\ {\lambda }_1^B\\ {\lambda }_{11}^{\mathrm{AB}}\end{array}\right],$(10) where the matrix with 0's and 1's is the design matrix formulating the loglinear model. Maximizing the likelihood gives maximum likelihood estimates ${\hat{\mu }}_{\mathrm{ij}}$ and parameter estimates $\hat{\lambda },{\hat{\lambda }}_1^A,{\hat{\lambda }}_1^B$ and ${\hat{\lambda }}_{11}^{\mathrm{AB}}$. In fitting the independence model, the last column of the design matrix, $(1,0,0,0{)}^T$ and the parameter ${\lambda }_{11}^{\mathrm{AB}}$ are left out of the equation. In software packages this can be accomplished in routines for Poisson regression, where the four observed counts $n_{\mathrm{ij}}$ are provided as the Y-variable, the model is represented by the design matrix in Equation (10(10) $\left[\begin{array}{c}\log ({\mu }_{11})\\ \log ({\mu }_{10})\\ \log ({\mu }_{01})\\ \log ({\mu }_{00})\end{array}\right]=\left[\begin{array}{cccc}1& 1& 1& 1\\ 1& 1& 0& 0\\ 1& 0& 1& 0\\ 1& 0& 0& 0\end{array}\right]\left[\begin{array}{c}\lambda \\ {\lambda }_1^A\\ {\lambda }_1^B\\ {\lambda }_{11}^{\mathrm{AB}}\end{array}\right],$(10) ), and the output consists of four fitted values ${\hat{\mu }}_{\mathrm{ij}}$ and the parameter estimates.

2.4.2. Generating two-way tables with given properties

We now develop a model for the two-way contingency table with prespecified marginal probabilities and a prespecified interaction parameter. For the interaction parameter we want to make use of the odds ratio ${\theta }^{\mathrm{AB}}$ defined in Equation (2(2) ${\theta }^{\mathrm{AB}}=\text{exp}({\lambda }_{11}^{\mathrm{AB}})=\frac{({\mu }_{11}/{\mu }_{01})}{({\mu }_{10}/{\mu }_{00})}=\frac{{\mu }_{11}{\mu }_{00}}{{\mu }_{01}{\mu }_{10}}.$(2) ). This can be accomplished by using an offset in Poisson regression. An offset is a column in the design matrix for which no parameter is estimated, i.e. the parameter is fixed to 1. The odds ratio is specified to give the desired level of dependence between lists. If we denote the prespecified odds ratio by ${\tilde{\theta }}^{\mathrm{AB}}$, then the desired loglinear interaction parameter is ${\tilde{\lambda }}_{11}^{\mathrm{AB}}=\log {\tilde{\theta }}^{\mathrm{AB}}$. Thus the population that we want to generate has to follow the following model: (11) $\log {\mu }_{\mathrm{ij}}=\lambda +{\lambda }_i^A+{\lambda }_j^B+{\tilde{\lambda }}_{\mathrm{ij}}^{\mathrm{AB}}.$(11) The following equation, an adjusted version of Equation (10(10) $\left[\begin{array}{c}\log ({\mu }_{11})\\ \log ({\mu }_{10})\\ \log ({\mu }_{01})\\ \log ({\mu }_{00})\end{array}\right]=\left[\begin{array}{cccc}1& 1& 1& 1\\ 1& 1& 0& 0\\ 1& 0& 1& 0\\ 1& 0& 0& 0\end{array}\right]\left[\begin{array}{c}\lambda \\ {\lambda }_1^A\\ {\lambda }_1^B\\ {\lambda }_{11}^{\mathrm{AB}}\end{array}\right],$(10) ), illustrates how this is implemented in Poisson regression: (12) $\left[\begin{array}{c}\log ({\mu }_{11})\\ \log ({\mu }_{10})\\ \log ({\mu }_{01})\\ \log ({\mu }_{00})\end{array}\right]=\left[\begin{array}{cccc}1& 1& 1& \log {\tilde{\theta }}^{\mathrm{AB}}\\ 1& 1& 0& 0\\ 1& 0& 1& 0\\ 1& 0& 0& 0\end{array}\right]\left[\begin{array}{c}\lambda \\ {\lambda }_1^A\\ {\lambda }_1^B\\ 1\end{array}\right].$(12) Notice that the parameter ${\lambda }_{11}^{\mathrm{AB}}$ from (10(10) $\left[\begin{array}{c}\log ({\mu }_{11})\\ \log ({\mu }_{10})\\ \log ({\mu }_{01})\\ \log ({\mu }_{00})\end{array}\right]=\left[\begin{array}{cccc}1& 1& 1& 1\\ 1& 1& 0& 0\\ 1& 0& 1& 0\\ 1& 0& 0& 0\end{array}\right]\left[\begin{array}{c}\lambda \\ {\lambda }_1^A\\ {\lambda }_1^B\\ {\lambda }_{11}^{\mathrm{AB}}\end{array}\right],$(10) ) is set to 1, i.e. it is not estimated, and this value 1 is multiplied with the prespecified $\log {\tilde{\theta }}^{\mathrm{AB}}$ in the design matrix. In software packages for Poisson regression this can be accomplished as follows:

• specify an offset vector which is $\log {\tilde{\theta }}^{\mathrm{AB}}$ when $(i,j)=(1,1)$ and 0 otherwise.

• use an independence model as the design matrix, as we only need estimates for $\lambda ,{\lambda }_i^A$ and ${\lambda }_j^B$.

• in order to fit a model that leads to estimated counts that follow the prespecified marginal probabilities, calculate joint probabilities for which these marginal probabilities hold, and use these joint probabilities, multiplied with some constant, as the counts $n_{\mathrm{ij}}$ to be analyzed with model (11(11) $\log {\mu }_{\mathrm{ij}}=\lambda +{\lambda }_i^A+{\lambda }_j^B+{\tilde{\lambda }}_{\mathrm{ij}}^{\mathrm{AB}}.$(11) ).

• then the estimates of expected frequencies can be used to derive the joint probabilities with the prespecified marginal probabilities and the prespecified odds ratio.

As an example, assume we want to generate population probabilities that have the prespecified marginal probabilities 0.7 and 0.8 and the prespecified odds ratio 2. Then the probabilities to be analyzed are 0.7 * 0.8 = 0.56, 0.7 * 0.2 = 0.14, 0.3 * 0.8 = 0.24 and 0.3 * 0.2 = 0.06, and multiply these probabilities with some number, such as 1000, to get counts. Thus 560, 140, 240 and 60 are the values of the dependent variable in Poisson regression. Specify $\log (2)$ as the offset and analyze the counts with the loglinear independence model and the offset. This yields the fitted values 584.8, 115.2, 215.2, 84.8 (rounded to 1 decimal place), and by dividing them by 1000 we have the four probabilities, with the prespecified marginal probabilities 0.7 and 0.8 and the prespecified odds ratio 2. The R-code we used can be found in the online Supplementary Material.

We now prove that for model (11(11) $\log {\mu }_{\mathrm{ij}}=\lambda +{\lambda }_i^A+{\lambda }_j^B+{\tilde{\lambda }}_{\mathrm{ij}}^{\mathrm{AB}}.$(11) ) the prespecified margins of the observed table are equal to those in the fitted table, irrespective of the offset. We follow Fienberg [27, Appendix II]. If we assume that the counts are observations from independent Poisson distributions, then the likelihood function is proportional to (13) $\prod _{\mathrm{ij}}{\mu }_{\mathrm{ij}}^{n_{\mathrm{ij}}}\exp (-{\mu }_{\mathrm{ij}}),$(13) and the kernel of the loglikelihood is, (14) $\sum _{\mathrm{ij}}{n}_{\mathrm{ij}}\log {\mu }_{\mathrm{ij}}-\sum _{\mathrm{ij}}{\mu }_{\mathrm{ij}}=\sum _{\mathrm{ij}}{n}_{\mathrm{ij}}(\lambda +{\lambda }_i^A+{\lambda }_j^B+{\tilde{\lambda }}_{\mathrm{ij}}^{\mathrm{AB}})-\sum _{\mathrm{ij}}{\mu }_{\mathrm{ij}}.$(14) At the maximum of the likelihood the derivative of the loglikelihood over the parameters is 0. So, for example, taking the derivative over ${\lambda }_i^A$ gives $n_{i+}-{\hat{\mu }}_{i+}=0$ and this gives $n_{i+}={\hat{\mu }}_{i+}$. Similarly for $n_{+j}={\hat{\mu }}_{+j}$. This completes the proof. It follows that including an offset for the interaction does not affect the property of ordinary loglinear models that for the fitted parameters the corresponding margins of the fitted values are equal to those of the observed values.

The model for the $3\times 3$ table in Section 2.1.1 has the following implementation: (15) $\left[\begin{array}{c}\log ({\mu }_{11})\\ \log ({\mu }_{12})\\ \log ({\mu }_{13})\\ \log ({\mu }_{21})\\ \log ({\mu }_{22})\\ \log ({\mu }_{23})\\ \log ({\mu }_{31})\\ \log ({\mu }_{32})\\ \log ({\mu }_{33})\end{array}\right]=\left[\begin{array}{cccccc}1& 1& 0& 1& 0& \log {\tilde{\theta }}_{11}^{\mathrm{AB}}\\ 1& 1& 0& 0& 1& \log {\tilde{\theta }}_{12}^{\mathrm{AB}}\\ 1& 1& 0& 0& 0& 0\\ 1& 0& 1& 1& 0& \log {\tilde{\theta }}_{21}^{\mathrm{AB}}\\ 1& 0& 1& 0& 1& \log {\tilde{\theta }}_{22}^{\mathrm{AB}}\\ 1& 0& 1& 0& 0& 0\\ 1& 0& 0& 1& 0& 0\\ 1& 0& 0& 0& 1& 0\\ 1& 0& 0& 0& 1& 0\end{array}\right]\left[\begin{array}{c}\lambda \\ {\lambda }_1^A\\ {\lambda }_2^A\\ {\lambda }_1^B\\ {\lambda }_2^B\\ 1\end{array}\right].$(15)

2.4.3. Model fitting for three-way tables

For the saturated model with three variables, the fitted values are equal to the observed counts so ${\hat{\mu }}_{\mathrm{ijk}}=n_{\mathrm{ijk}}$. As an example of a restricted model, the estimates for the [AB][AC] model can be calculated directly as (16) ${\hat{\mu }}_{\mathrm{ijk}}=\frac{{\hat{\mu }}_{\mathrm{ij}+}{\hat{\mu }}_{i+k}}{{\hat{\mu }}_{i++}}=\frac{n_{\mathrm{ij}+}{n}_{i+k}}{n_{i++}}.$(16) However, for the loglinear model [AB][AC][BC] estimates must be calculated indirectly, iteratively (for example, [6], p. 84).

We now develop (4(4) $\log {\mu }_{\mathrm{ijk}}=\lambda +{\lambda }_i^A+{\lambda }_j^B+{\lambda }_k^C+{\lambda }_{\mathrm{ij}}^{\mathrm{AB}}+{\lambda }_{\mathrm{ik}}^{\mathrm{AC}}+{\lambda }_{\mathrm{jk}}^{\mathrm{BC}}+{\lambda }_{\mathrm{ijk}}^{\mathrm{ABC}},$(4) ) in the form of Poisson regression. There are eight elements $n_{\mathrm{ijk}}$, being counts from eight Poisson distributions. The likelihood can be set up for ${\mu }_{\mathrm{ijk}}$ and maximised over the parameters. Model (4(4) $\log {\mu }_{\mathrm{ijk}}=\lambda +{\lambda }_i^A+{\lambda }_j^B+{\lambda }_k^C+{\lambda }_{\mathrm{ij}}^{\mathrm{AB}}+{\lambda }_{\mathrm{ik}}^{\mathrm{AC}}+{\lambda }_{\mathrm{jk}}^{\mathrm{BC}}+{\lambda }_{\mathrm{ijk}}^{\mathrm{ABC}},$(4) ) can be denoted as, (17) $\left[\begin{array}{c}\log ({\mu }_{111})\\ \log ({\mu }_{110})\\ \log ({\mu }_{101})\\ \log ({\mu }_{100})\\ \log ({\mu }_{011})\\ \log ({\mu }_{010})\\ \log ({\mu }_{001})\\ \log ({\mu }_{000})\end{array}\right]=\left[\begin{array}{cccccccc}1& 1& 1& 1& 1& 1& 1& 1\\ 1& 1& 1& 0& 1& 0& 0& 0\\ 1& 1& 0& 1& 0& 1& 0& 0\\ 1& 1& 0& 0& 0& 0& 0& 0\\ 1& 0& 1& 1& 0& 0& 1& 0\\ 1& 0& 1& 0& 0& 0& 0& 0\\ 1& 0& 0& 1& 0& 0& 0& 0\\ 1& 0& 0& 0& 0& 0& 0& 0\end{array}\right]\left[\begin{array}{c}\lambda \\ {\lambda }_1^A\\ {\lambda }_1^B\\ {\lambda }_1^C\\ {\lambda }_{11}^{\mathrm{AB}}\\ {\lambda }_{11}^{\mathrm{AC}}\\ {\lambda }_{11}^{\mathrm{BC}}\\ {\lambda }_{111}^{\mathrm{ABC}}\end{array}\right].$(17) As before, for (17(17) $\left[\begin{array}{c}\log ({\mu }_{111})\\ \log ({\mu }_{110})\\ \log ({\mu }_{101})\\ \log ({\mu }_{100})\\ \log ({\mu }_{011})\\ \log ({\mu }_{010})\\ \log ({\mu }_{001})\\ \log ({\mu }_{000})\end{array}\right]=\left[\begin{array}{cccccccc}1& 1& 1& 1& 1& 1& 1& 1\\ 1& 1& 1& 0& 1& 0& 0& 0\\ 1& 1& 0& 1& 0& 1& 0& 0\\ 1& 1& 0& 0& 0& 0& 0& 0\\ 1& 0& 1& 1& 0& 0& 1& 0\\ 1& 0& 1& 0& 0& 0& 0& 0\\ 1& 0& 0& 1& 0& 0& 0& 0\\ 1& 0& 0& 0& 0& 0& 0& 0\end{array}\right]\left[\begin{array}{c}\lambda \\ {\lambda }_1^A\\ {\lambda }_1^B\\ {\lambda }_1^C\\ {\lambda }_{11}^{\mathrm{AB}}\\ {\lambda }_{11}^{\mathrm{AC}}\\ {\lambda }_{11}^{\mathrm{BC}}\\ {\lambda }_{111}^{\mathrm{ABC}}\end{array}\right].$(17) ) the matrix with 0's and 1's is the design matrix formulating the loglinear model. A restricted model can be fitted by dropping columns in Equation (17(17) $\left[\begin{array}{c}\log ({\mu }_{111})\\ \log ({\mu }_{110})\\ \log ({\mu }_{101})\\ \log ({\mu }_{100})\\ \log ({\mu }_{011})\\ \log ({\mu }_{010})\\ \log ({\mu }_{001})\\ \log ({\mu }_{000})\end{array}\right]=\left[\begin{array}{cccccccc}1& 1& 1& 1& 1& 1& 1& 1\\ 1& 1& 1& 0& 1& 0& 0& 0\\ 1& 1& 0& 1& 0& 1& 0& 0\\ 1& 1& 0& 0& 0& 0& 0& 0\\ 1& 0& 1& 1& 0& 0& 1& 0\\ 1& 0& 1& 0& 0& 0& 0& 0\\ 1& 0& 0& 1& 0& 0& 0& 0\\ 1& 0& 0& 0& 0& 0& 0& 0\end{array}\right]\left[\begin{array}{c}\lambda \\ {\lambda }_1^A\\ {\lambda }_1^B\\ {\lambda }_1^C\\ {\lambda }_{11}^{\mathrm{AB}}\\ {\lambda }_{11}^{\mathrm{AC}}\\ {\lambda }_{11}^{\mathrm{BC}}\\ {\lambda }_{111}^{\mathrm{ABC}}\end{array}\right].$(17) ). For example, in fitting the model where B and C are independent given A, the last two columns of the design matrix and the last two parameters are left out of the equation.

2.4.4. Generating three-way tables with given properties

We now focus on loglinear models with prescribed odds ratios. We start with the model with no three-factor interaction. The model for a three-way contingency table with three prespecified odds ratios ${\tilde{\theta }}^{\mathrm{AB}},{\tilde{\theta }}^{\mathrm{AC}}$ and ${\tilde{\theta }}^{\mathrm{BC}}$ is denoted as, (18) $\log {\mu }_{\mathrm{ij}}=\lambda +{\lambda }_i^A+{\lambda }_j^B+{\lambda }_k^C+\log {\tilde{\theta }}^{\mathrm{AB}}+\log {\tilde{\theta }}^{\mathrm{AC}}+\log {\tilde{\theta }}^{\mathrm{BC}}.$(18) In the adjusted version of (17(17) $\left[\begin{array}{c}\log ({\mu }_{111})\\ \log ({\mu }_{110})\\ \log ({\mu }_{101})\\ \log ({\mu }_{100})\\ \log ({\mu }_{011})\\ \log ({\mu }_{010})\\ \log ({\mu }_{001})\\ \log ({\mu }_{000})\end{array}\right]=\left[\begin{array}{cccccccc}1& 1& 1& 1& 1& 1& 1& 1\\ 1& 1& 1& 0& 1& 0& 0& 0\\ 1& 1& 0& 1& 0& 1& 0& 0\\ 1& 1& 0& 0& 0& 0& 0& 0\\ 1& 0& 1& 1& 0& 0& 1& 0\\ 1& 0& 1& 0& 0& 0& 0& 0\\ 1& 0& 0& 1& 0& 0& 0& 0\\ 1& 0& 0& 0& 0& 0& 0& 0\end{array}\right]\left[\begin{array}{c}\lambda \\ {\lambda }_1^A\\ {\lambda }_1^B\\ {\lambda }_1^C\\ {\lambda }_{11}^{\mathrm{AB}}\\ {\lambda }_{11}^{\mathrm{AC}}\\ {\lambda }_{11}^{\mathrm{BC}}\\ {\lambda }_{111}^{\mathrm{ABC}}\end{array}\right].$(17) ) this becomes (using that $\log {\tilde{\theta }}^{\mathrm{AB}}+\log {\tilde{\theta }}^{\mathrm{AC}}+\log {\tilde{\theta }}^{\mathrm{BC}}=\log {\tilde{\theta }}^{\mathrm{AB}}{\tilde{\theta }}^{\mathrm{AC}}{\tilde{\theta }}^{\mathrm{BC}}$): (19) $\left[\begin{array}{c}\log ({\mu }_{111})\\ \log ({\mu }_{110})\\ \log ({\mu }_{101})\\ \log ({\mu }_{100})\\ \log ({\mu }_{011})\\ \log ({\mu }_{010})\\ \log ({\mu }_{001})\\ \log ({\mu }_{000})\end{array}\right]=\left[\begin{array}{ccccc}1& 1& 1& 1& \log {\tilde{\theta }}^{\mathrm{AB}}{\tilde{\theta }}^{\mathrm{AC}}{\tilde{\theta }}^{\mathrm{BC}}\\ 1& 1& 1& 0& \log {\tilde{\theta }}^{\mathrm{AB}}\\ 1& 1& 0& 1& \log {\tilde{\theta }}^{\mathrm{AC}}\\ 1& 1& 0& 0& 0\\ 1& 0& 1& 1& \log {\tilde{\theta }}^{\mathrm{BC}}\\ 1& 0& 1& 0& 0\\ 1& 0& 0& 1& 0\\ 1& 0& 0& 0& 0\end{array}\right]\left[\begin{array}{c}\lambda \\ {\lambda }_1^A\\ {\lambda }_1^B\\ {\lambda }_1^C\\ 1\end{array}\right],$(19) where the parameter for the offsets has been set to 1. In software packages for Poisson regression this can be accomplished as follows:

• specify $\log {\tilde{\theta }}^{\mathrm{AB}}$, $\log {\tilde{\theta }}^{\mathrm{AC}}$ and $\log {\tilde{\theta }}^{\mathrm{BC}}$, and use them to generate the offset column of the design matrix as in (19(19) $\left[\begin{array}{c}\log ({\mu }_{111})\\ \log ({\mu }_{110})\\ \log ({\mu }_{101})\\ \log ({\mu }_{100})\\ \log ({\mu }_{011})\\ \log ({\mu }_{010})\\ \log ({\mu }_{001})\\ \log ({\mu }_{000})\end{array}\right]=\left[\begin{array}{ccccc}1& 1& 1& 1& \log {\tilde{\theta }}^{\mathrm{AB}}{\tilde{\theta }}^{\mathrm{AC}}{\tilde{\theta }}^{\mathrm{BC}}\\ 1& 1& 1& 0& \log {\tilde{\theta }}^{\mathrm{AB}}\\ 1& 1& 0& 1& \log {\tilde{\theta }}^{\mathrm{AC}}\\ 1& 1& 0& 0& 0\\ 1& 0& 1& 1& \log {\tilde{\theta }}^{\mathrm{BC}}\\ 1& 0& 1& 0& 0\\ 1& 0& 0& 1& 0\\ 1& 0& 0& 0& 0\end{array}\right]\left[\begin{array}{c}\lambda \\ {\lambda }_1^A\\ {\lambda }_1^B\\ {\lambda }_1^C\\ 1\end{array}\right],$(19) ).

• use an independence model in the design matrix, as we only need estimates for $\lambda ,{\lambda }_i^A$, ${\lambda }_j^B$ and ${\lambda }_k^C$.

• in order to fit a model that leads to estimated counts that follow prespecified marginal probabilities, calculate joint probabilities for which these marginal probabilities hold by multiplying the relevant probabilities. Use these joint probabilities, multiplied by some constant representing a notional population size, as the observed counts corresponding with $n_{\mathrm{ijk}}$ to be analyzed with model (19(19) $\left[\begin{array}{c}\log ({\mu }_{111})\\ \log ({\mu }_{110})\\ \log ({\mu }_{101})\\ \log ({\mu }_{100})\\ \log ({\mu }_{011})\\ \log ({\mu }_{010})\\ \log ({\mu }_{001})\\ \log ({\mu }_{000})\end{array}\right]=\left[\begin{array}{ccccc}1& 1& 1& 1& \log {\tilde{\theta }}^{\mathrm{AB}}{\tilde{\theta }}^{\mathrm{AC}}{\tilde{\theta }}^{\mathrm{BC}}\\ 1& 1& 1& 0& \log {\tilde{\theta }}^{\mathrm{AB}}\\ 1& 1& 0& 1& \log {\tilde{\theta }}^{\mathrm{AC}}\\ 1& 1& 0& 0& 0\\ 1& 0& 1& 1& \log {\tilde{\theta }}^{\mathrm{BC}}\\ 1& 0& 1& 0& 0\\ 1& 0& 0& 1& 0\\ 1& 0& 0& 0& 0\end{array}\right]\left[\begin{array}{c}\lambda \\ {\lambda }_1^A\\ {\lambda }_1^B\\ {\lambda }_1^C\\ 1\end{array}\right],$(19) ).

• then, as before in Section 2.4.1, the estimates of the expected frequencies can be used to derive the joint probabilities with the prespecified marginal probabilities and the prespecified odds ratios.

As an example, assume we want to generate population probabilities that have the prespecified marginal probabilities 0.8, 0.7 and 0.9 for variables A, B and C respectively, and with all the prespecified odds ratios $=2$. The probabilities following the independence model are 0.7 * 0.8 * 0.9 = 0.504, 0.7 * 0.8 * 0.1 = 0.056, 0.7 * 0.2 * 0.9 = 0.216, 0.7 * 0.2 * 0.1 = 0.024, 0.3 * 0.8 * 0.9 = 0.126, 0.3 * 0.8 * 0.1 = 0.014, 0.3 * 0.2 * 0.9 = 0.054 and 0.3 * 0.2 * 0.1 = 0.006 and we multiply these probabilities by 1,000 to get counts 504, 56, 216, 24, 126, 14, 54 and 6 which form the values of the dependent variable in Poisson regression. Specify $\log (2)$ as the three offsets and analyze the counts with the loglinear independence model and the offset as in Equation (19(19) $\left[\begin{array}{c}\log ({\mu }_{111})\\ \log ({\mu }_{110})\\ \log ({\mu }_{101})\\ \log ({\mu }_{100})\\ \log ({\mu }_{011})\\ \log ({\mu }_{010})\\ \log ({\mu }_{001})\\ \log ({\mu }_{000})\end{array}\right]=\left[\begin{array}{ccccc}1& 1& 1& 1& \log {\tilde{\theta }}^{\mathrm{AB}}{\tilde{\theta }}^{\mathrm{AC}}{\tilde{\theta }}^{\mathrm{BC}}\\ 1& 1& 1& 0& \log {\tilde{\theta }}^{\mathrm{AB}}\\ 1& 1& 0& 1& \log {\tilde{\theta }}^{\mathrm{AC}}\\ 1& 1& 0& 0& 0\\ 1& 0& 1& 1& \log {\tilde{\theta }}^{\mathrm{BC}}\\ 1& 0& 1& 0& 0\\ 1& 0& 0& 1& 0\\ 1& 0& 0& 0& 0\end{array}\right]\left[\begin{array}{c}\lambda \\ {\lambda }_1^A\\ {\lambda }_1^B\\ {\lambda }_1^C\\ 1\end{array}\right],$(19) ). This yields the fitted values (rounded to 1 decimal place) 547.3, 39.5, 186.3, 26.9, 99.0, 14.3, 67.4 and 19.4. By dividing these by 1000 we have the eight probabilities with the prespecified marginal probabilities 0.7, 0.8 and 0.9 and the three prespecified odds ratios 2. For example, the marginal probability of being in variable A equals the sum of the fitted values for A = 1 divided by 1000, which is 0.8 and the conditional odds ${\theta }^{\mathrm{BC}|A=1}=(547.3\ast 26.9)/(39.5\ast 186.3)=2$. The R-code we used can be found in the online Supplementary Material.

We now discuss the model with a prespecified three-factor interaction. Above we discussed model fitting for three-way tables, with no three-factor interaction. Here we develop the model to include a prespecified three-factor interaction parameter in a loglinear model. The three-factor interaction is defined as the state where two conditional odds ratios are not equal, for example, ${\tilde{\theta }}^{\mathrm{AB}|C=0}\ne {\tilde{\theta }}^{\mathrm{AB}|C=1}$. These conditional odds ratios are related to loglinear parameters by ${\tilde{\theta }}^{\mathrm{AB}|C=0}=\exp {\tilde{\lambda }}_{11}^{\mathrm{AB}}$ and ${\tilde{\theta }}^{\mathrm{AB}|C=1}=\exp [{\tilde{\lambda }}_{11}^{\mathrm{AB}}+{\tilde{\lambda }}_{111}^{\mathrm{ABC}}]$. As similar results hold for ${\tilde{\theta }}^{\mathrm{AC}|B=0}$, ${\tilde{\theta }}^{\mathrm{AC}|B=1}$, ${\tilde{\theta }}^{\mathrm{BC}|A=0}$ and ${\tilde{\theta }}^{\mathrm{BC}|A=1}$ this shows that the three-way interaction is represented by $\exp {\tilde{\lambda }}_{111}^{\mathrm{ABC}}$. A three-way interaction will be handled as follows. In Equation (19(19) $\left[\begin{array}{c}\log ({\mu }_{111})\\ \log ({\mu }_{110})\\ \log ({\mu }_{101})\\ \log ({\mu }_{100})\\ \log ({\mu }_{011})\\ \log ({\mu }_{010})\\ \log ({\mu }_{001})\\ \log ({\mu }_{000})\end{array}\right]=\left[\begin{array}{ccccc}1& 1& 1& 1& \log {\tilde{\theta }}^{\mathrm{AB}}{\tilde{\theta }}^{\mathrm{AC}}{\tilde{\theta }}^{\mathrm{BC}}\\ 1& 1& 1& 0& \log {\tilde{\theta }}^{\mathrm{AB}}\\ 1& 1& 0& 1& \log {\tilde{\theta }}^{\mathrm{AC}}\\ 1& 1& 0& 0& 0\\ 1& 0& 1& 1& \log {\tilde{\theta }}^{\mathrm{BC}}\\ 1& 0& 1& 0& 0\\ 1& 0& 0& 1& 0\\ 1& 0& 0& 0& 0\end{array}\right]\left[\begin{array}{c}\lambda \\ {\lambda }_1^A\\ {\lambda }_1^B\\ {\lambda }_1^C\\ 1\end{array}\right],$(19) ) we fill in $\log {\tilde{\theta }}^{\mathrm{AB}|C=0}$ in cell (2,5) of the design matrix and $\log {\tilde{\theta }}^{\mathrm{AB}|C=1}$ in cell (1,5) of the design matrix, and similarly for the other conditional odds ratios. Thus, the model with prespecified conditional odds ratios can be presented as, (20) $\left[\begin{array}{c}\log ({\mu }_{111})\\ \log ({\mu }_{110})\\ \log ({\mu }_{101})\\ \log ({\mu }_{100})\\ \log ({\mu }_{011})\\ \log ({\mu }_{010})\\ \log ({\mu }_{001})\\ \log ({\mu }_{000})\end{array}\right]=\left[\begin{array}{ccccc}1& 1& 1& 1& \log {\tilde{\theta }}^{\mathrm{AB}|C=1}{\tilde{\theta }}^{\mathrm{AC}|B=1}{\tilde{\theta }}^{\mathrm{BC}|A=1}\\ 1& 1& 1& 0& \log {\tilde{\theta }}^{\mathrm{AB}|C=0}\\ 1& 1& 0& 1& \log {\tilde{\theta }}^{\mathrm{AC}|B=0}\\ 1& 1& 0& 0& 0\\ 1& 0& 1& 1& \log {\tilde{\theta }}^{\mathrm{BC}|A=0}\\ 1& 0& 1& 0& 0\\ 1& 0& 0& 1& 0\\ 1& 0& 0& 0& 0\end{array}\right]\left[\begin{array}{c}\lambda \\ {\lambda }_1^A\\ {\lambda }_1^B\\ {\lambda }_1^C\\ 1\end{array}\right].$(20) We can trivially write (21) ${\tilde{\theta }}^{\mathrm{AB}|C=1}={\tilde{\theta }}^{\mathrm{AB}|C=0}\frac{{\tilde{\theta }}^{\mathrm{AB}|C=1}}{{\tilde{\theta }}^{\mathrm{AB}|C=0}},$(21) the term $[{\tilde{\theta }}^{\mathrm{AB}|C=1}/{\tilde{\theta }}^{\mathrm{AB}|C=0}]$, which is the difference in odds for C = 1 and C = 0, is equal to $\exp {\tilde{\lambda }}_{111}^{\mathrm{ABC}}$. The terms ${\tilde{\theta }}^{\mathrm{AC}|B=1}$ and ${\tilde{\theta }}^{\mathrm{BC}|A=1}$ include this same three-way interaction term, i.e. (22) $\begin{array}{rl}{\tilde{\theta }}^{\mathrm{AC}|B=1}& ={\tilde{\theta }}^{\mathrm{AC}|B=0}\frac{{\tilde{\theta }}^{\mathrm{AB}|C=1}}{{\tilde{\theta }}^{\mathrm{AB}|C=0}},\end{array}$(22) (23) $\begin{array}{rl}{\tilde{\theta }}^{\mathrm{BC}|A=1}& ={\tilde{\theta }}^{\mathrm{BC}|A=0}\frac{{\tilde{\theta }}^{\mathrm{AB}|C=1}}{{\tilde{\theta }}^{\mathrm{AB}|C=0}}.\end{array}$(23) In fact since the three-factor interaction parameter ${\tilde{\lambda }}_{111}^{\mathrm{ABC}}$ is the difference between each pair of conditional odds ratios, it is clear that (24) $\frac{{\theta }^{\mathrm{AB}|C=1}}{{\theta }^{\mathrm{AB}|C=0}}=\frac{{\theta }^{\mathrm{AC}|B=1}}{{\theta }^{\mathrm{AC}|B=0}}=\frac{{\theta }^{\mathrm{BC}|A=1}}{{\theta }^{\mathrm{BC}|A=0}}.$(24) If we specify the denominators of each fraction, then as soon as one numerator is specified then all the fractions are determined. There are therefore only four free odds ratios in Equation (24(24) $\frac{{\theta }^{\mathrm{AB}|C=1}}{{\theta }^{\mathrm{AB}|C=0}}=\frac{{\theta }^{\mathrm{AC}|B=1}}{{\theta }^{\mathrm{AC}|B=0}}=\frac{{\theta }^{\mathrm{BC}|A=1}}{{\theta }^{\mathrm{BC}|A=0}}.$(24) ), and the three-factor interaction can be specified by providing all the denominators and any one of the numerators. Therefore we can rewrite element (1,5) of the design matrix in (20(20) $\left[\begin{array}{c}\log ({\mu }_{111})\\ \log ({\mu }_{110})\\ \log ({\mu }_{101})\\ \log ({\mu }_{100})\\ \log ({\mu }_{011})\\ \log ({\mu }_{010})\\ \log ({\mu }_{001})\\ \log ({\mu }_{000})\end{array}\right]=\left[\begin{array}{ccccc}1& 1& 1& 1& \log {\tilde{\theta }}^{\mathrm{AB}|C=1}{\tilde{\theta }}^{\mathrm{AC}|B=1}{\tilde{\theta }}^{\mathrm{BC}|A=1}\\ 1& 1& 1& 0& \log {\tilde{\theta }}^{\mathrm{AB}|C=0}\\ 1& 1& 0& 1& \log {\tilde{\theta }}^{\mathrm{AC}|B=0}\\ 1& 1& 0& 0& 0\\ 1& 0& 1& 1& \log {\tilde{\theta }}^{\mathrm{BC}|A=0}\\ 1& 0& 1& 0& 0\\ 1& 0& 0& 1& 0\\ 1& 0& 0& 0& 0\end{array}\right]\left[\begin{array}{c}\lambda \\ {\lambda }_1^A\\ {\lambda }_1^B\\ {\lambda }_1^C\\ 1\end{array}\right].$(20) ) as (25) $\log {\tilde{\theta }}^{\mathrm{AB}|C=1}{\tilde{\theta }}^{\mathrm{AC}|B=1}{\tilde{\theta }}^{\mathrm{BC}|A=1}=\log {\tilde{\theta }}^{\mathrm{AB}|C=0}{\tilde{\theta }}^{\mathrm{AC}|B=0}{\tilde{\theta }}^{\mathrm{BC}|A=0}{\left[\frac{{\tilde{\theta }}^{\mathrm{AB}|C=1}}{{\tilde{\theta }}^{\mathrm{AB}|C=0}}\right]}^3,$(25) using only four odds ratios.

As an example, assume that we want to create population probabilities with conditional odds ratios ${\tilde{\theta }}^{\mathrm{AB}|C=0}=2$ and ${\tilde{\theta }}^{\mathrm{AB}|C=1}=3$. It follows that the three factor term is $[{\tilde{\theta }}^{\mathrm{AB}|C=1}/{\tilde{\theta }}^{\mathrm{AB}|C=0}]=3/2$. Then, let ${\tilde{\theta }}^{\mathrm{AC}|B=1}=1$, and it follows that ${\tilde{\theta }}^{\mathrm{AC}|B=0}[{\tilde{\theta }}^{\mathrm{AB}|C=1}/{\tilde{\theta }}^{\mathrm{AB}|C=0}]=1$, so that ${\tilde{\theta }}^{\mathrm{AC}|B=0}=2/3$. Similarly, if ${\tilde{\theta }}^{\mathrm{BC}|A=1}=1$, then ${\tilde{\theta }}^{\mathrm{BC}|A=0}=2/3$. Thus the model equation becomes (26) $\left[\begin{array}{c}\log ({\mu }_{111})\\ \log ({\mu }_{110})\\ \log ({\mu }_{101})\\ \log ({\mu }_{100})\\ \log ({\mu }_{011})\\ \log ({\mu }_{010})\\ \log ({\mu }_{001})\\ \log ({\mu }_{000})\end{array}\right]=\left[\begin{array}{ccccc}1& 1& 1& 1& \log 3\\ 1& 1& 1& 0& \log 2\\ 1& 1& 0& 1& \log (2/3)\\ 1& 1& 0& 0& 0\\ 1& 0& 1& 1& \log (2/3)\\ 1& 0& 1& 0& 0\\ 1& 0& 0& 1& 0\\ 1& 0& 0& 0& 0\end{array}\right]\left[\begin{array}{c}\lambda \\ {\lambda }_1^A\\ {\lambda }_1^B\\ {\lambda }_1^C\\ 1\end{array}\right].$(26)

3. Application to population size estimation

3.1. Dual system estimation

Nowadays dual system estimation methods are commonly used to estimate population totals for domains of interest. The dual system estimator is the simplest capture-recapture model and is used to estimate the part of the population that is missing after having linked individuals in two overlapping lists. Originally this method was implemented to estimate the size and density of wildlife populations, but it has also been developed and used in official statistics, where the dual-system estimator was originally used for the 1950 decennial United States Census to estimate coverage error [28].

The linked data can be presented in a two-way contingency table: we observe the counts of those in both lists, those in list A only, and those in list B only, but not those missing from both lists. Thus there are three observed cell counts, and the missing count for those in neither of the lists is to be estimated, which can conveniently be done using the loglinear model framework [29].

Dual system estimation relies on five key assumptions, outlined by the International Working Group for Disease Monitoring and Forecasting [30]:

1. The population is closed, i.e. there is no change such as births or deaths.

2. There is perfect matching between the two lists.

3. The capture/inclusion probability of elements in at least one of the lists is homogeneous [31,32].

4. The probabilities of inclusion in the lists are independent.

5. There are no erroneous captures in any list (no overcoverage).

One of the key assumptions is that the probabilities of inclusion in the lists are independent, which is easily violated [24]. However, with the data used for population size estimation, this assumption is impossible to verify and so the design of both lists is very important, to create the conditions where this assumption is likely to hold. A way to reduce the impact of violating this strict assumption is to divide the population by certain characteristics, such as age, sex and some geographical breakdown, so that we need only assume that the lists are independent conditional on these classifiers.

Assumption 4 follows on from assumption 3, which implies that individuals who where included in the first list and those who were not have the same probability of being included in the second list, so that presence in the first list does not affect response in the second list, and the lists are therefore statistically independent [30]. Therefore, as there are three observed counts, to estimate cell counts ${\hat{\mu }}_{\mathrm{ij}}$ for a two-way contingency table only three parameters can be used and no interaction term can be modelled, and (1(1) $\log {\mu }_{\mathrm{ij}}=\lambda +{\lambda }_i^A+{\lambda }_j^B+{\lambda }_{\mathrm{ij}}^{\mathrm{AB}},$(1) ) is adjusted to include no interaction term, i.e. ${\lambda }_{\mathrm{ij}}^{\mathrm{AB}}=0$. Independence implies that the odds ratio (2(2) ${\theta }^{\mathrm{AB}}=\text{exp}({\lambda }_{11}^{\mathrm{AB}})=\frac{({\mu }_{11}/{\mu }_{01})}{({\mu }_{10}/{\mu }_{00})}=\frac{{\mu }_{11}{\mu }_{00}}{{\mu }_{01}{\mu }_{10}}.$(2) ) is 1, so that the missing count is estimated by ${\hat{\mu }}_{00}={\hat{\mu }}_{10}{\hat{\mu }}_{01}/{\hat{\mu }}_{11}=n_{10}{n}_{01}/n_{11}$.

3.2. Triple system estimation

The triple system estimator uses three linked lists. Counts can be summarized in a three-way contingency table where one count is missing, namely the count for the individuals that are missed by all three lists. By estimating this count and adding it to the seven observed counts, the population size is estimated. Unlike a two-way contingency table, in a three-way contingency table the pairwise interactions between lists can be modelled. This means that the homogeneous capture/inclusion probabilities and the independence assumption can be relaxed, while the other assumptions outlined for dual system estimation still hold.

As there are seven cells with observed counts, only loglinear models with up to seven independent parameters can be fitted. This is usually solved by assuming the three-factor interaction term ${\hat{\lambda }}_{\mathrm{ijk}}^{\mathrm{ABC}}=0$. Including two-factor interaction terms in the model enables us to estimate population totals where pairwise relationships are present [6]. For this paper, we are interested in fitting the maximal model for a three-way contingency table with one unobserved cell, which is an adjusted version of (4(4) $\log {\mu }_{\mathrm{ijk}}=\lambda +{\lambda }_i^A+{\lambda }_j^B+{\lambda }_k^C+{\lambda }_{\mathrm{ij}}^{\mathrm{AB}}+{\lambda }_{\mathrm{ik}}^{\mathrm{AC}}+{\lambda }_{\mathrm{jk}}^{\mathrm{BC}}+{\lambda }_{\mathrm{ijk}}^{\mathrm{ABC}},$(4) ) with no three factor interaction term.

3.3. How to apply the proposed methodology to population size estimation

The following outlines how to generate contingency tables with prespecified margins and interactions between lists. Such generated contingency tables can be used in simulations to study the behaviour of estimators for the population size, as the true population size is known.

To generate a contingency table where lists are independent, the total population size and inclusion probabilities for each list are prespecified. The cell counts are then computed by multiplying the marginal inclusion probabilities and the population size. This results in cell counts that follow statistical independence between lists.

To generate a table with a given dependence structure, an independent loglinear model is then fitted to these cell counts with the inclusion of an offset term, as presented in (11(11) $\log {\mu }_{\mathrm{ij}}=\lambda +{\lambda }_i^A+{\lambda }_j^B+{\tilde{\lambda }}_{\mathrm{ij}}^{\mathrm{AB}}.$(11) ) and (18(18) $\log {\mu }_{\mathrm{ij}}=\lambda +{\lambda }_i^A+{\lambda }_j^B+{\lambda }_k^C+\log {\tilde{\theta }}^{\mathrm{AB}}+\log {\tilde{\theta }}^{\mathrm{AC}}+\log {\tilde{\theta }}^{\mathrm{BC}}.$(18) ). Then the fitted counts have the following properties:

1. Their sum is equal to the prespecified population size,

2. The prespecified odds ratios hold, and

3. The prespecified marginal response probabilities hold.

In order to investigate the behaviour of population size estimators, multinomial samples are drawn from this population. For the study of the dual system estimator, in each sample the cell (0,0) is set to missing and the estimator is applied. For the study of the triple systems estimator, in each sample the cell (0,0,0) is set to missing and the estimator is applied. When t samples are drawn, t population size estimates are obtained that can be compared with the prespecified population size.

3.4. Simulation study

The aim of the simulation study is to illustrate the methodology. Probabilities are generated using prespecified marginal probabilities and odds ratios. 2000 multinomial samples are drawn with these probabilities and with the total population size set equal to 1000 or 2000. The prespecified probabilities and odds ratios mirror those used by Brown et al. [22] in their investigation of methods to deal with dependence in the coverage estimation for the 2001 population census of England and Wales. The outcomes of the simulations are presented in Tables 1 and 2. The results include the specified values for the total population, the marginal response probabilities, the odds ratio (OR) between the lists, the mean and median of the population size estimates, empirical 95% confidence intervals (which is the range between the 2.5% and 97.5% quantiles of the estimates), the relative bias %, $(\hat{\mu }-N)/N$ and the Coefficient of Variation $\hat{\sigma }$/N, where N is the prespecified population size, $\hat{\mu }$ is the mean and $\hat{\sigma }$ is the standard deviation of the population size estimates over the simulations.

Table 1. Simulation results for dual systems estimation with inclusion probabilities ${\pi }_A$ and ${\pi }_B$, sizes of odds ratios (OR), the mean and median of the estimates, the 95% confidence interval of the mean, the relative bias and coefficient of variation (CV), where the fitted model is [A][B].

N	${\pi }_A$	${\pi }_B$	OR	Mean	Median	95%CI		rbias%	CV
1000	0.8	0.7	1	1000.3	1000.0	979.5,	1020.6	0.027	0.011
1000	0.9	0.8	1	1000.0	1000.0	989.8,	1010.0	−0.005	0.005
2000	0.8	0.7	1	2000.1	2000.2	1971.3,	2029.1	0.005	0.007
2000	0.9	0.8	1	2000.1	2000.1	1985.2,	2014.6	0.007	0.004
1000	0.8	0.7	1.1	994.5	994.5	973.9,	1015.2	−0.545	0.011
1000	0.9	0.8	1.1	998.0	998.1	986.9,	1008.4	−0.197	0.005
2000	0.8	0.7	1.1	1989.9	1989.1	1960.3,	2018.0	−0.553	0.007
2000	0.9	0.8	1.1	1996.1	1996.0	1980.9,	2010.2	−0.197	0.004
1000	0.8	0.7	2	957.8	957.8	937.5,	978.9	−4.215	0.011
1000	0.9	0.8	2	984.4	984.6	972.3,	995.9	−1.560	0.006
2000	0.8	0.7	2	1915.6	1915.7	1886.4,	1943.8	−4.219	0.008
2000	0.9	0.8	2	1968.6	1968.7	1950.8,	1985.5	−1.571	0.004

Table 2. Simulation results for triple systems estimation with inclusion probabilities ${\pi }_A$ = 0.8, ${\pi }_B$ = 0.7 and ${\pi }_C$ = 0.9, with different sizes of odds ratios (OR), specified models, the mean and median of the estimates, the 95% confidence interval of the mean, the relative bias and coefficient of variation (CV).

N	OR (AB,AC,BC)	Model	Mean	Median	95%CI		rbias%	CV
1000	(1, 1, 1)	[A][B][C]	1000.0	1000.2	994.7	1004.6	0.002	0.003
1000	(1, 1, 1)	[AB][AC][BC]	1000.2	1000.1	993.5	1007.6	0.020	0.004
2000	(1, 1, 1)	[A][B][C]	2000.1	2000.1	1992.9	2007.1	0.003	0.002
2000	(1, 1, 1)	[AB][AC][BC]	2000.2	2000.1	1991.3	2010.6	0.012	0.002
1000	(1.5, 2, 1)	[A][B][C]	993.7	993.9	986.7	1000.1	−0.627	0.003
1000	(1.5, 2, 1)	[AB][AC][BC]	1000.4	1000.0	990.0	1013.0	0.036	0.006
1000	(1.5, 2, 1)	[AB][AC]	1000.2	1000.0	991.2	1009.4	0.018	0.005
1000	(1.5, 2, 1)	[A][BC]	993.0	993.2	986.0	999.4	−0.704	0.004
2000	(1.5, 2, 1)	[A][B][C]	1987.4	1987.4	1977.7	1996.9	−0.629	0.002
2000	(1.5, 2, 1)	[AB][AC][BC]	2000.2	1999.6	1985.8	2017.9	0.009	0.004
2000	(1.5, 2, 1)	[AB][AC]	2000.1	2000.2	1987.6	2014.0	0.007	0.003
2000	(1.5, 2, 1)	[A][BC]	1985.9	1985.9	1976.2	1995.5	−0.706	0.003

3.4.1. Two-way table

We set the response probabilities to either, for list A (${\pi }_A$) 0.8 and for list B (${\pi }_B$) 0.7; or for list A (${\pi }_A$) 0.9 and for list B (${\pi }_B$) 0.8. The joint probabilities (${\pi }_{11},{\pi }_{10},{\pi }_{01},{\pi }_{00}$) are, respectively, (0.56, 0.24, 0.14, 0.06) and (0.72, 0.18, 0.08 and 0.02). Cell counts are obtained by ${\mu }_{\mathrm{ij}}=(N{\pi }_{\mathrm{ij}})$.

The results presented in Table 1 from these simulations show, that when the lists are independent (OR = 1), and the independence model is fitted, the mean of the estimated population size across 2000 multinomial samples is very close to the true population size. However, the bias appears to depend on the size of the population and the response probabilities for each list. When the total population is smaller, N = 1000, and the response probabilities are lower, ${\pi }_A$ = 0.7 and ${\pi }_B$ = 0.8, for the lists, the mean of the estimated population sizes is slightly larger than the true population size. However, the bias is negligible when N = 1000 and ${\pi }_A$ = 0.8 and ${\pi }_B$ = 0.9, as the mean estimated population total across all multinomial samples is approximately equal to the true population total. Therefore, the higher the response probabilities of the lists, the less biased the estimates of the population size will be.

When the lists are dependent, and the independence model is fitted, the mean estimated population size is smaller than the true population size. This shows that positive odds ratios result in population size estimates being negatively biased (compare with [24]). In some cases, when the odds ratio = 1.1, the true population size falls within the confidence interval. This is an issue of power, as the misspecification of the model is small, and the number of simulations is apparently not large enough to detect it.

3.4.2. Three-way table

In our simulation, the three lists have response probabilities 0.8 for list A (${\pi }_A$), 0.7 for list B (${\pi }_B$) and 0.9 for list C (${\pi }_C$). The joint probabilities (${\pi }_{111},{\pi }_{110},{\pi }_{101},{\pi }_{100}$, ${\pi }_{011},{\pi }_{010},{\pi }_{001},{\pi }_{000}$) are, respectively, (0.504, 0.056, 0.216, 0.024, 0.126, 0.054, 0.014, 0.006). Cell counts are obtained by ${\mu }_{\mathrm{ijk}}=(N{\pi }_{\mathrm{ijk}})$.

The results of the simulation are shown in Table 2. The model used to generate the population data is specified by the odds ratios in the column OR, and the model used to analyze the data in the column Model. Thus there are correctly specified models, such as OR = (1,1,1) and Model is [A][B][C], overparameterised models, such as OR = (1,1,1) and Model is [AB][AC][BC], and misspecified models, such as OR = (1.5, 2, 1) and Model is [A][BC]. When the model is specified correctly the relative bias is close to zero. The variance of the estimates, decreases as the population size increases, compare, for example, OR = (1,1,1) and Model [A][B][C] with population sizes 1000 and 2000. When the model is overparameterised, the variance is larger than when the model is correctly specified.

4. Discussion

We proposed a simple method to generate contingency tables with fixed marginal probabilities and dependence structures described by loglinear models. We discuss the generation of two-way and three-way contingency tables and this can be extended to in a straightforward way.

We specify the loglinear model in terms of Poisson regression with a log link. Using an offset, prespecified interaction parameters can be fitted that can be understood in terms of odds ratios. For two-way tables, odds ratios describe dependence between two variables, and conditional odds ratios describe partial dependence between three or more variables in higher-way tables.

We think that this method is simpler than previously proposed methods which are discussed in the introduction, and it can be extended to other models that can be presented as loglinear models, such as logistic regression models with categorical predictors and some latent variable models. The structure of our approach to generating data is the same as the standard approach for undertaking multiple system estimation, so it is straightforward to understand how the dependence is introduced. Our approach is less tailored to the generalized estimating equations which have been the main application of other approaches.

In this paper we apply this method to population size estimation for two-way and three-way contingency tables. The results show the impacts of dependence between lists and model (mis)specification.

Disclosure statement

No potential conflict of interest was reported by the author(s).