Abstract

11
Introduction
In epidemiological studies, observed data are often collected subject to misclassification errors. Such misclassification errors cause bias of estimation and reduce efficiency in the analysis of contingency tables. Many investigators have discussed how to adjust for the effects of misclassi- fication, see the review papers, Chen (1989) and Walter and Irwig (1987).  Espeland and Odoroff (1985) discussed the maximum likelihood estimation for a recursive system of log-linear models based on double sampling schemes by the EM algorithm.  Espeland and Hui (1987) applied the Fisher scoring algorithm to evaluating variances and covariances of estimates. From the Bayesian viewpoint, Geng and Asano (1989) considered estimation methods for misclassified categorical data making use of prior information and double sampling schemes, and obtained posterior means as estimates.  Viana (1994) applied Bayesian computations based on the matrix of misclassification probabilities to small-sample multinomial data. Evans, Guttman, Hatiovsky and Swartz (1996) discussed the implementation of the Gauss-Jacobi quadrature and the Gibbs sampling algorithm for the posterior analysis of binary response data with misclassification. In this paper, we present a Bayesian approach that utilizes prior knowledge about misclas- sification and incorporates this prior knowledge with observations subject to misclassification. Although the EM algorithm or the Fisher scoring algorithm are often applied to estimating model parameters, these algorithms can not evaluate posterior distributions on the model pa- rameters.  Furthermore, these algorithms do not apply to our estimation problem because of the unidentifiability of the model parameters.  However, our Bayesian approach, assuming a prior distribution on the model parameters, can overcome these problems. In order to find the posterior distribution of model parameters and calculate posterior means as estimates of them, we use the data augmentation(DA) algorithm by Tanner and Wong (1987). In Section 2, we show the Bayesian computation to find a posterior distribution given mis- classified observed data. In Section 3, we give the DA algorithm to approximate the posterior distribution, because of difficulties with the calculation of the posterior distribution. Section 4 presents two numerical experiments to examine the performance of the DA algorithm.
2
Misclassified observed data and Bayesian computation
Let X and Y be categorical variables having I and J categories, respectively, and let Y be a misclassified variable of Y having K categories. We assume that two types of misclassified data are observed: (i) data for X and Y', denoted as n = {ni+k | i ∈ {1,….,I},k ∈ {1,….,K}}, and (ii) data for Y and Y ′ , denoted as m = { m + jk | j ∈ {1, . . , J },k ∈ {1, . . , K }}, where the symbol "+" means the sum over corresponding variables, for example, n i+k = ∑ j n ijk. Let p ijk denote a probability for (X,Y,Y,Y') = (i,j,k) and θχγγ, = {P ijk | i ∈ {1,….,I},j ∈ {1, . . , J}, k ∈ {1, . . , K}} denote a set of probabilities. In this paper, the goal is to find the posterior distribution of model parameters θχγ = { Pij + | i ∈ {1,...,I},j ∈ {1,...,J}} which are the marginal probabilities of X and Y , and obtain the posterior means of θ χ  as estimates. Assume that n and m have independently multinomial distributions with parameters θ γγ = {Pi+k | i ∈ {1,...,I},k ∈ {1,...,K}}, and θγγ , = {P+jk | j ∈ {1,...,J},k ∈ {1,...,K}}, respectively, that is,
and that the prior distribution of θ χγγ  is a Dirichlet distribution which has the density function
where αχγγ, = {αijk | i ∈ {1,….,I},j ∈ {1,….,J},k ∈ {1,….,K}} is a set of hyper- parameters of the prior distribution of θχγγι.  From the equations (1) and (2), we obtain the mixture posterior distribution given n and m . The posterior density is given by
where ∑ Ω(n) denotes the sum over all possible { ñijk } under the conditions ñijk ≥ 0 for all i, j and k , and ∑ j ñ ijk = π i+k , and ∑ Ω(m) denotes the sum over all possible { ñ ijk } under the conditions ˜h ijk ≥ 0 for all i, j and k , and ∑ i ˜ujk = m +jk . However, the posterior density (3) has a very complicated function, because of combinational explosion. It is extremely difficult to calculate exactly the posterior distribution and these calculation may take a long time when observed data are moderately large. For the case that the incomplete-data posterior density is complicated as the equation (3) and the complete-data posterior density is relative easy to handle and draw from, the DA algorithm is very suitable. In next section, we present the DA algorithm to approximate the posterior distribution to find the posterior distribution and estimate posterior means of model parameters θχγ,
3
DA algorithm for misclassified data
The DA algorithm consists of iterating between the imputation-step and the posterior-step. For this misclassified multinomial model, the DA algorithm is given by the following iterative scheme: Imputation-step:  Repeat the following two steps for l = 1, . . , L to o obtain the imputed data of n and m such as ñ = {ñijk | i ∈ {1,...,I},j ∈ {1,...,J},k ∈ {1,...,I}, ∑j ñijk = πi+k} and ñ = { ñ ijk | i ∈ {1, . . ,I}, j ∈ {1, . . ,J}, k ∈ {1, . . ,K}, ∑ i ñ ijk = m +jk}. 1. Generate cell probabilities { p ijk } from the current estimate of the posterior distribu- tion, 2. Generate the imputed data ñ(l) and ñ(l) from the predictive distributions which have the conditional multinomial distributions, given n and m , with densities
where P j | i,k = P ijk/P i+k and P i | j,k = P ijk/P +jk Posterior-step:  Update the current approximation of the posterior distribution of θχγγ, given these imputed data ñ(l) and ñ(l), for l = 1,..., L, by the Monte Carlo method,
Until the approximated distribution converges to a stationary distribution, the imputation-step and the posterior-step are iterated.  Getting a stationary distribution, the values of L may be increased to improve the accuracy with respect to the Monte Carlo method. Then we can find the posterior distribution of θχγ and can obtain easily the posterior means and variances of them. Furthermore, it is possible to calculate the highest posterior density(HPD) region that is the Bayesian analog of the confidence intervals.
4
Numerical experiments
We provide two numerical experiments to examine the performance of the DA algorithm de- scribed the previous section.
Comparison of estimates and exact values
4.1
In the following numerical experiment, we compare the estimates obtained by the DA algo- rithm with the posterior means using the exact Bayesian calculation given by Geng and Asano (1987) regarding the data from Diamond and Lilienfeld (1962) that reported a case-control s- tudy concerning the circumcision status of male partners of woman with and without cervical cancer. The study sample was categorized by cervical cancer status, X (Case and Control), and self-reported circumcision status, Y ′ (Yes or No), in the left side of Table 1. In order to gain the information on the degree of misclassification of circumcision status, the supplemen- tal sample concerning the relationship between actual circumcision status, Y (Yes or No), and Y  was gathered from the separate population shown by in the center of Table 1.  Espeland and Hui (1987) described that, for the misclassified multinomial model, the conditional inde- pendence model between X and Y given Y ′ from the class of hierarchical log-linear models was appropriate, since no observed data for X , Y  and Y ′ were obtained. Furthermore, for the conditional independence model, Geng and Asano (1989) gave the hypothetical prior informa- tion shown in the right side of Table 1.  We use their prior information as hyper-parameters αχγγ, = {α j k | i ∈ {Case, Control}, j ∈ {Yes, No}, k ∈ {Yes, No}} and then obtain a posterior distribution and estimates of θχγ = {pij+ | i ∈ {Case, Control},j ∈ {Yes, No}} for X and Y by posterior means. In this numerical experiment, we e aluate the accuracy of the estimates using the DA algo- rithm in comparison with the exact posterior means given by Geng and Asano (1989). Table 2 shows the exact values, and the posterior means, the standard deviations(SDs) and the posterior 95 % credible intervals(CIs) of θ χγ obtained from simulated 10,000 samples after a burn-in of 1,000 samples. It can be seen that the estimates have approximately three-digit accuracy for the exact values. From the numerical results, we can see that the DA algorithm works quite well to estimate posterior means.
Performance of the DA algorithm
4.2
In the next numerical experiment, we e xamine the performance of the DA algorithm in compari- son with the EM algorithm and the Fisher scoring algorithm. We apply the DA algorithm to the double sampling data from Hochberg (1977). The data were the highway safety research data relating the seat-belt usages to driver injuries. The main sample was of 80,084 accidents that were recorded by police subject to misclassification errors. The subsample was of 1,796 accidents that were recorded by both imprecise police reports and precise hospital interviews. Then, by the double sampling design, the subsample was randomly selected from the main sample. Thus, the subsample and the main sample have independent and identical distributions. The main sample and the subsample in Table 3 were categorized by four variables X, X ′ , Y and Y ′ , where X and Y denote precise personal survey for seat-belt usages and driver injuries, and X ′ and Y ′ denote imprecise police reports for them. In this experiment, we estimate model parameters under the saturated multinomial model, because our purpose is to investigate whether the DA algorithm is applicable to estimate model parameters, but not to analyze the misclassified observed data. For these data, we assume that the main sample data and the subsample data have inde- pendent and identical multinomial distributions with θ XייYY = {Pijkl | i ∈ {Yes, No}, j ∈ {Yes, No}, k ∈ {Yes, No}, l ∈ {Yes, No}}, where p ijkl = Pr(X = i,X' = j,Y = k,Y' = l) and the prior distribution for θχיציhas the Dirichlet distribution with hyper-parameters αχχιγγι = {αιjkι | i ∈ {Yes, No},j ∈ {Yes, No}, k ∈ {Yes, No}, l ∈ {Yes, No}}. Then the model parameters are marginal probabili-
where p i+k+ = ∑ j,t P ijk . We utilize the subsample in Table 3 as hyper-parameters α χ ιγγ , and obtain estimates of θχγ by the DA algorithm. Table 4 shows the estimates and the SDs of θχγ obtained by the DA algorithm, the exact Bayesian calculation, the Fisher scoring algorithm and the EM algorithm.  The estimates using the DA algorithm can be found from simulated samples 100,000 after a burn-in samples 10,000 in two chains.  The exact values of estimates of θχγ using the Bayesian calculation are given by Geng and Asano (1989) who assumed the Jeffreys noninformative prior. The estimation using the Fisher scoring algorithm were carried out with ℓ EM developed by Vermunt (1997). From these numerical results, it can be seen that the DA algorithm has the equivalent perfor- mance of the EM and the Fisher scoring algorithm in comparison with these estimates and SDs. Applying the EM and the Fisher scoring algorithm, their algorithms have such disadvantages as it is impossible to find the posterior distribution of model parameters, may not be applied them to estimation of parameters owing to unidentifiability of the model and may be difficult to calculate the Fisher information matrix needed in the Fisher scoring algorithm.
Concluding remarks
55
In this paper, we discussed the DA algorithm to estimate model parameters for misclassified categorical data. We gave the posterior distribution by exact Bayesian computation. To avoid complicated calculation, we used the DA algorithm and find the posterior distribution.  It is easily seen that the DA algorithm is the iterative simulation version of the EM algorithm that the imputation-step corresponds to the E-step and the posterior-step corresponds to the M-step. In order to explore the possibility of parameter estimation by the DA algorithm, we provided two numerical experiments.  In the first experiment, we evaluated accuracy of estimates in comparison with exact values.  In the second experiment, we e examined the performance of the DA algorithm. The results of both the two numerical experiments showed the advantage of applying the DA algorithm in terms of accuracy of estimates and in terms of algorithm simplicity to find the posterior distribution. For the inference of multidimensional contingency tables, the Bayesian inference by the DA algorithm can be e asily extended and also widely utilized.  Then we may need to take account for conditional independence between variables in models. For parameters assuming the conditional independence model, A prior Dirichlet distribution that has hyper Markov laws by Dawid and Lauritzen (1993) is very suitable.  Our future problem is how to incorporate prior information with the hyper-parameters of a hyper Dirichlet prior distribution without consistency.
Acknowledgements
The authors would like to thank an editor and a referee for helpful suggestions and correcting our paper.