Subjective Elicitation of Dirichlet Hyperparameters Using Past Data : A Study of Ovarian Cancer Patients

Elicitation of prior plays a very important role in Bayesian paradigm especially when dealing with rare disease problems in medical field. The reason being that we do not get enough data to draw valid inferences always. Since the subject of study is human population, one cannot do experiments with their health. The prior distribution supports the final results by some additional information gained from the experts. In any case if an appropriate expert is not available, we can use past data to get information about the prior and its hyperparameters. The present paper provides a technique of elicitation of prior hyperparameters based on a well known multinomialDirichlet model. Since the main focus is on medical data problems, the inferences on odds ratios and interaction parameters are also provided. Numerical illustration is based on a real dataset from Israel on patients having ovarian cancer. Although the details have been given in the context of ovarian cancer patients, the development in the paper is equally well applicable for any such disease.


Introduction
Bayesian approaches to the analysis of epidemiological data represent a powerful tool for interpretation of study results and evaluation of hypotheses about exposure-disease relations.Definitely prior distributions have a major role in Bayesian inference and one has to be very careful while choosing them.A small error in the choice of prior may result into misleading inferences.There are various types of priors suggested in the literature, say, for example, Jefferys' prior, reference prior, conjugate prior, uniform prior, etc.In true sense, the spirit of Bayesian paradigm lies certainly in the subjective elicitation of prior that may be approached by the experimenter through her subjective assessment or expert opinion.The objective is to convert the assessment or the opinion into a probabilistic form.In most of the cases the expert is unaware of statistical terminologies.As a result, it becomes a difficult task to extract the information required to draw inferences on a particular hypothesis.Moreover, an elicitation is considered appropriate if the distribution that is derived accurately represents the expert's knowledge regardless of how good that knowledge is.The idea of subjective prior elicitation has been advocated by a number of researchers (see, for example, O'Hagan (2006)) but not given due emphasis by most of the applied Bayesian practitioners.
Elicitation or subjective elicitation generally involves two strategies.In one case it may involve elic-iting the form of the prior distribution whereas in the other case it may be concerned with eliciting only the hyperparameters of an assumed form of prior.The situation may also be considered where one requires the elicitation of both prior distribution and its hyperparameters.This paper focusses on elicitation of hyperparameters assuming the multinomial-Dirichlet combination considered earlier by several authors.This may be equivalently considered as specifying the exact member of the Dirichlet family based on the expert opinion that may be appropriate for the considered situation.Since the prior distribution does play a crucial role in final inferences, a small ambiguity in expert's perception or a small misinterpretation of subjective opinion in converting into an appropriate a priori judgement may result into drastically poor inferences.Keeping this in mind, a number of strategies have been suggested in the literature for eliciting the hyperparameters of a chosen prior distribution, in general, and of a conjugate prior distribution, in particular.Staël von Holstein (1971) is perhaps the first reference where the author suggested a method for assessing the conjugate prior for a Bernoulli process.His technique for assessing a beta distribution made use of the estimates of median and first and third quartiles from the experts' subjective opinion.Kadane, Dickey, Winkler, Smith, and Peters (1980) presented a method for estimating conjugate priors for a linear regression model.Their method made use of the multivariate t predictive distribution, which involved the assessment of belief about measures of central tendency for the regression coefficients as well as an assessment of belief about variation or covariation and the appropriate value of a degrees-of-freedom parameter.Besides, we also have an important reference by Winkler, Smith, and Kulkarni (1978) where authors considered the elicitation of conjugate priors for the linear regression model by using predictive distributions.
Among other significant references, Chaloner and Duncan (1983) presented a method called predictive modal estimation in the context of assessing beta prior.In their method, the expert first provides an assessment of the mode of the distribution and then assesses the likelihood of other points along the distribution relative to the likelihood of the mode.Similarly, Garthwaite and Dickey (1988) examined the linear regression problem but introduced a technique that was based on the concept of points of constrained minimum variance.In this technique, certain values of the independent variables are given to the expert.The expert is then asked to select values for the remaining independent variables so that his or her uncertainty regarding the dependent variable is minimized.A few fractile assessments are also required from the expert in this approach to complete the elicitation of the subjective probability distributions.
Truly speaking, to determine the prior hyperparameters, expert judgement is usually sought on the quantities that may be easily assessed.The elicited quantities are then equated to their theoretical expressions in order to solve for the unknown prior hyperparameters.A few such quantities may include the central tendency measures such as the mean, median, mode, etc. Quantile estimates may also be sometimes used based on the expert judgement.The problem, however, encountered with most of these elicited characteristics is that they result in numerically challenging scenarios when equated with the corresponding theoretical expressions.Dorp and Mazzuchi (2003) is an important reference that makes use of the quantile estimates for specification of the parameters of the beta distribution and its multivariate analogue.They have given a detailed accountability of the inherent numerical intricacies and accordingly provided a few possible solutions based on certain iterative procedures.This paper considers a problem from genetic epidemiology where it is supposed that disease in a human population occurs sometimes due to one's physical structure inherited by birth and sometimes due to outer environmental components that become part of individual's everyday life.Subtle differences in genetic factors also cause people to respond differently to the same environmental exposure.The complex interaction between an individual's genetic composition and environmental agents is the cause for almost all the diseases.Although the interaction of genes and environment are often discussed in the context of disease or negative traits, the impact of their interaction can be protective, neutral, or even harmful.Mukherjee and Chatterjee (2008), Mukherjee, Ahn, Gruber, Ghosh, and Chatterjee (2010) and Gupta and Upadhyay (2014) are some recent references in this area where the authors have considered geneenvironment association problems and obtained inferences on odds ratios and various interaction parameters.The current work provides a complete Bayes analysis of the same problem considering a multinomial-Dirichlet modelling assumption.It mainly focusses on the elicitation of Dirichlet prior hyperparameters by the method given by Dorp and Mazzuchi (2003).As we do not have expert knowledge on the specific application discussed in this paper, we presume that this information is available from other sources.It is important to mention here that the objective of the present paper is to formalize the procedure presuming as if we have the expert opinion though in actual implementation we never got an expert who could be examined for her opinion.As such, we rely on an alternative strategy based on the past data or on a small subset of the given data to derive the expert opinion on the quartiles that are subsequently used in the elicitation of prior hyperparameters.There is no claim made here that our procedure is perfect.We believe, however, that our development can serve as the beginning of a new research where prior is given due consideration in the light of expert's opinion to get an appropriate conclusion.
The plan of the paper is as follows.Section 2 discusses the formulation of the model for the problem under consideration with a focus on our main objective, the subjective elicitation of prior.The complete method is described in subsection 2.1.Based on the method discussed in this subsection, the numerical illustration is provided in Section 3.Here we consider an ovarian cancer dataset on Israeli women and accordingly draw the relevant inferences.We also provide the same study on reduced dataset and the results are presented in subsection 3.1.Finally, a brief conclusion is given at the end.

Model formulation
At the outset let us consider a specific epidemiological problem that involves case-control scenario in the structure of 2 × 2 × (l + 1).This structure is relevant where there are exactly three categorizations with each of the first two categorizations, F 1 and F 2 , having two levels and the third categorization, F 3 , having (l + 1) levels.We may refer, for instance, the two levels of F 2 , that is, F 2 = 0 and F 2 = 1 as controls and cases, respectively.The complete structure consisting of n(= n 0 + n 1 ) observations can be classified as in Table 1 in the form of different cell counts.We have used the notation r ijk to represent the count of individual cell where both i and j correspond to F 1 and F 2 , respectively, and k corresponds to the third categorization, F 3 .Obviously, both i and j are binary variables taking values either 0 or 1 whereas the variable k is a polytomous variable that may take values 0, 1, ..., l.
It may be noted that we have considered two multinomial distributions, one corresponding to controls (F 2 = 0) and the other corresponding to cases (F 2 = 1).Moreover, the components of p 0 and p 1 are the corresponding cell probabilities given by p ijk = r ijk /n j where Σr ijk = n j , i = 0, 1; j = 0, 1; k = 0, 1, ..., l.The likelihood functions based on the above two configurations for controls and cases can be written as Let us now consider separate Dirichlet priors for the cell probabilities p j = (p 0j0 , ..., p 0jl , p 1j0 , ..., p 1jl ) for j = 0, 1 as where a j = (a 0j0 , ..., a 0jl , a 1j0 ..., a 1jl ), j = 0, 1, are the hyperparameters.The corresponding posteriors up to proportionality can be written as where the posterior for the controls (F 2 = 0) corresponds to j = 0 and the posterior for the cases (F 2 = 1) corresponds to j = 1.Obviously, the two posteriors given in (3), when normalized, are again Dirichlet distributions with updated parameters (r j + a j ) where r j and a j for j = 0, 1 are already defined above.
In epidemiological studies, the objective generally lies in drawing inferences about the odds ratio, a ratio that measures the odds of exposure for cases against controls.These odds ratios can be calculated for any of the factors of interest, say, F 1 , F 2 , F 3 and also for the joint effect of these factors, say, for example, effect of F 2 and F 3 on different levels of F 1 , etc. Let OR F 3k = p 000 p 01k /p 00k p 010 denotes the odds ratio associated with the k th level of F 3 for F 1 = 0 and OR F 1 = p 000 p 110 /p 010 p 100 denotes the odds ratio associated with F 1 for F 3 = 0.There may be situations when any two of these three components might have a joint affect on the third one.To check whether such association exists, we define the measure of association between, say, F 1 and F 3 considering only its k th level and keeping where the subscript k with F 3 denotes its k th level.If, however, this association comes out to be non-zero, one may go a step ahead for calculating the multiplicative interaction parameter between F 1 and F 3 at its k th level (see also Mukherjee and Chatterjee (2008)).This interaction parameter can be given as ψ k = (p 00k p 010 p 100 p 11k )/(p 000 p 01k p 10k p 110 ), k = 1, ..., l. (5)

Subjective elicitation of prior
To present the methodology in a simple way, we consider the same modelling formulation given earlier.We, however, use slightly different form in this subsection and the symbols used here may not be exactly related to what we have defined earlier.To begin with, let us consider x = (x 1 , . . ., x ), p = (p 1 , . . ., p ), λ = (λ 1 , . . ., λ ) and define To describe briefly the elicitation procedure for the model under consideration, let us first reparameterize the Dirichlet (λ) distribution by introducing new parameters β = i λ i and α i = λ i /β, i = 1, ..., , yielding the density, Thus to obtain the estimates of λ i s, we need to elicit both α i s and β.
It is well known that beta distribution is a special case of Dirichlet distribution and all the marginal distributions of Dirichlet variables follow beta distributions.An important property of beta distributions is that for all p ∈ [0, 1], 0 < m 1 < m 2 ⇒ g(p|m 1 , n) > g(p|m 2 , n) and n 1 > n 2 > 0 ⇒ g(p|m, n 1 ) > g(p|m, n 2 ) where p denotes a univariate analogue of p and g(p|m, n) denotes a beta distribution with parameters m and n given by If we re-parameterize the beta distribution in the same way as it was suggested for the Dirichlet distribution, the above property can be written as, where P is used to denote the random variable corresponding to p and α (.) and β are the parameters of the re-parameterized beta distribution.Obviously, the marginals of Dirichlet distribution also follow this property.Finally, the quantile constraint concept given below in the form of a definition can be used as well (see, for example, Dorp and Mazzuchi (2003)).
This property can be easily extended for the Dirichlet distribution at the level of each of its marginals.Thus to get the values of α i , i = 1, 2, ..., , and β, we need to solve α and β for p ∼ g(p|α, β) given in ( 6) under the two quantile constraints (p i qL , q i L ) and (p i qU , q i U ) for any P i , where q i L < q i U , and ( − 1) single quantile constraint (p j q , q j ) for P j , j = 1, 2, ..., , j = i.As usual, P (.) is used to denote the random variable corresponding to the component p (.) .Since this problem cannot be solved in a closed form, one has to resort to an appropriate numerical solution.A numerical solution is defined below (see also Dorp and Mazzuchi (2003)).A word of remark: although the method of Dorp and Mazzuchi (2003) is straightforward to deal with, it suffers from an important drawback in the sense that P i 's are not treated symmetrically.The method actually requires placing two quantile constraints on one of the probabilities and ( −1) single quantile constraint on others.Moreover, P is treated quite differently than other probabilities (see, for example, Evans, Guttman, and Li (2017)).Alternative methods are, of course, suggested in the literature (see Zapata-Vázquez, O'Hagan, and Bastos (2014) and Evans et al. (2017)) to deal with the elicitation of Dirichlet parameters efficiently but we stick to the method of Dorp and Mazzuchi (2003) simply because of its inherent ease.Moreover, since our illustration in the next section considers a large data size, it is expected that a slight deviation in elicited prior will not affect our final inferences.
To clarify the details, let us suppose an initial interval (c 1 , f 1 ) that contains β * (an estimate of β) and let an initial value of β * be β * 1 = c 1 +f 1 2 .For this particular β * 1 , calculate α * i1 , an initial value of α i , by considering an interval (d 1 , e 1 ) = (0, 1), which is supposed to contain α * i1 and take , where p qU is the upper quantile value suggested by the expert.If (q U ) 1 ≤ q U , where q U is the area of upper quartile, then (d 2 , e 2 ) ≡ (d 1 , α * i1 ), otherwise (d 2 , e 2 ) ≡ (α * i1 , e 1 ).This follows from the property of the beta distribution given by (8).Now take α * i2 = d 2 +e 2 2 and again calculate (q U ) 2 = P r(P ≤ p qU |α * i2 , β * 1 ) and repeat the same procedure described above until at the t th stage (q U ) t ≈ q U .After this we shall go on for refining β * 1 .For this, we calculate q L = P r(P ≤ (p qL ) 1 |α * it , β * 1 ), where q L is the area of lower quantile.If (p qL ) 1 < p qL , where p qL is the lower quantile value given by the expert, then For this value β * 2 of β again update the value of α i by the same procedure.The value of β is updated until (p qL ) t ≈ p qL .In this way the final values of α i and β are obtained as α * i and β * , respectively.Using the final value of β * , we calculate α * j , j = 1, 2, ..., , j = i through the same procedure by which α * i is calculated.To determine the initial interval (c 1 , f 1 ) containing β * , first set c 1 = 0 as β > 0. To obtain the upper bound f 1 of this initial interval, set β 11 = 1 which implies f 1 = 2. Now solve for α i1 satisfying (q U ) 1 = P r(P ≤ p qU |α i1 , β 11 ).Proceeding in a similar manner described above, we may go on for updating α i1 .Say, for α it at t th stage, we get (q U ) t ≈ q U .We may then find out (p qL ) 1,t by solving q L = P r(P ≤ (p qL ) 1,t |α it , β 1,t ).In case (p qL ) 1,t < p qL then set β 1,t+1 = 2 * β 1,t .We may repeat the above procedure for all t.Conversely, if (p qL ) 1,t > p qL , set f 1 = β 1,t .In this way the initial interval containing β * may be determined.

Numerical illustration
To illustrate our procedure, we considered an ovarian cancer dataset based on the women of Israel.This is a partially real and a partially simulated dataset which has been analyzed earlier by a number of authors (see, for example, Modan, Hartge, Hirsh-Yechezkel, Chetrit, Lubin, Beller, Ben-Baruch, Fishman, Menczer, Struewing, Tucker, and Wacholder (2001), Chatterjee and Carroll (2005), Mukherjee and Chatterjee (2008), etc.).A recent analysis of this data has been done by Gupta and Upadhyay (2014) but in an empirical Bayes way.The dataset is given in Table 2.The categories F 1 , F 2 and F 3 are denoted by G, D and E representing genetic susceptibility (absent/ present), disease status (controls/ cases) and environmental exposure (absent/ present at level k, k = 1, 2, 3), respectively.The gene responsible for the occurrence of ovarian cancer is supposed to be BRCA1 and/or 2 mutation and the levels of environmental exposures, considered as oral contraceptive (OC) use and parity, are defined as follows: • For OC use: E = 0 ⇒ subjects who never used OC, E = 1 ⇒ corresponds to those who used OC up to 3 years, E = 2 ⇒ corresponds to those who used OC from 3 years to 6 years and E = 3 ⇒ corresponds to those who used OC for more than 6 years.
• For parity: E = 0 ⇒ corresponds to women who have no children, E = 1 ⇒ corresponds to women having 1-3 children, E = 2 ⇒ corresponds to women having 3-6 children and E = 3 ⇒ corresponds to women having more than 6 children.
For convenience most of the previous authors have analyzed this problem considering that genes and environmental components are independent of each other.This dataset can be converted into a 2 × 2 × 2 data by combining the non-zero cells of the category environmental exposure to E = 1.Such an attempt might be of general interest to those medical practitioners who simply want to study the impact of presence or absence of environmental components and are not interested to go into a detailed study at different levels.
We also considered a dataset of moderate sample size with n = 50 extracted from 2 × 2 × 2 setup discussed above in such a way that each cell frequency remains at least equal to unity.The dataset is shown in Table 3 and it is taken exclusively for the purpose of comparison.It is to be noted that such datasets with small to very small sample sizes may not result in the desired inferences on association and interaction parameters as the corresponding reported cell frequencies may not be the true representatives of the prevalence of gene and environmental components in the actual population.
In our analysis the primary objective is to estimate the odds ratio, measure of association between G and E and, in case this association is non-zero, the multiplicative interaction parameter between the latter two.It may be noted, however, that the G − E association in the control population at k th level of E (see ( 4)) and the multiplicative interaction parameter between G and E = k, k = 1, 2, 3, (see ( 5)) based on the observations in Table 2 are not logically appealing as a number of cell frequencies are too small to draw any valid conclusion on the actual G − E association or the multiplicative interaction parameter.We, therefore, propose to consider such measures for 2 × 2 × 2 setup along with 2 × 2 × 4 setup.In this case, we shall drop the subscript k from OR E k , θ GE k , and ψ k keeping their interpretations the same.
Since we did not have any expert who can be contacted to get her subjective opinion, we divided the entire data into two parts, one having first 50 observations and other, the remaining 1529 observations.This division was done for both the 2 × 2 × 4 and 2 × 2 × 2 setups separately.The first part of the data was considered as past data and it was used for elicitation of hyperparameters whereas the second part of the data was used for the desired inferences.The two datasets are shown in Tables 4-5.The values in the parentheses correspond to 2 × 2 × 2 setup.
Table 4: Past data with sample size 50 for elicitation of prior hyperparameters Table 5: Remaining data with sample size 1529 for the inferential developments  In order to elicit the prior hyperparameters (see subsection 2.1), we first obtained the empirical Bayes (EB) estimates of Dirichlet hyperparameters a j , j = 0, 1, where a 0 denotes the vector of hyperparameters corresponding to controls and a 1 denotes the same for cases, based on the past data of size 50 (see also Gupta and Upadhyay (2014)).It is to be noted that some of the cell observations were zero in past data making the corresponding cell probabilities also zero.In order to find the EB estimates for these cases, we considered very small cell probabilities of the order 10 −5 for these cells so that the corresponding EB estimates can be worked out.Before we proceed further, let us describe briefly the EB procedure used in the present paper.
The EB approach traditionally uses a prior density with hyperparameters estimated on the basis of observed data usually by means of classical approach.If maximum likelihood (ML) estimates are easily obtainable, the same are often preferred.To provide a brief outline of used EB procedure, let us consider x = (x 1 , . . ., x ) as the set of observed multinomial data.The corresponding multinomial parameters (p 1 , . . ., p ) can be easily obtained as the estimates of different cell probabilities based on various cell counts.Let these estimates are denoted by p = (p 1 , . . ., p ) (see also Gupta and Upadhyay (2014)).Finally, the EB estimates of the parameters of the Dirichlet distribution, which is the prior for multinomial parameters, can be obtained by maximizing the corresponding log-likelihood function given by, There are a number of methods for numerically maximizing this objective function G(λ) as there is no closed form solution for the same.A detailed survey for various methods can be found in Lwin and Maritz (1989), Minka (2000), etc. (see also Gupta and Upadhyay (2014)).We, however, used the procedure given by Minka (2000).The results of EB estimates are shown in Table 6 where the bracketed values represent the EB estimates corresponding to 2 × 2 × 2 setup.It may be noted here that we have described the EB procedure based on the model formulation given at the beginning of subsection 2.1 although the procedure was actually implemented on (2) for j = 0, 1.The description given on ( 9) is simply meant for notational convenience.
Table 6: EB estimates of the components of a j (j = 0, 1) using past data No doubt, the EB criterion is close to classical paradigm because of the involvement of observed data in its evaluation.The method is often criticized by subjective Bayesians who consider prior information exogenous to observations.We, therefore, finally resort to subjective elicitation of hyperparameters based on the description given in subsection 2.1.Since the method is based on elicitation of hyperparameters using quartiles obtained from the experts, we actually work for these quartiles using EB estimates derived from the past data due to non-availability of experts in our case.To be explicit, because we do not have an expert available to specify the quartiles used in the elicitation, we estimated the hyperparameters of the Dirichlet distribution using the EB method just described and applied to a subset of the data.This Dirichlet gives the quartiles that are subsequently used in the elicitation procedure described in subsection 2.1.We do not claim that relying on EB estimates is the only possibility.One can, of course, use alternative method as well.Our next step, therefore, consisted of evaluating the quartiles by solving the incomplete beta functions (see subsection 2.1) after replacing the unknown Dirichlet parameters by their EB estimates given in Table 6 (see also (8) and the Definition given in subsection 2.1).Our final step consisted of estimating the parameters α i , i = 1, ..., 8(4) and β by the method suggested in subsection 2.1.These estimated parameters are shown separately for OC use and parity in Table 7 and 8 where Table 7 corresponds to 2 × 2 × 4 setup and Table 8 corresponds to the same for 2 × 2 × 2 setup.These estimates then provide the estimates of original Dirichlet hyperparameters a j , j = 0, 1, used in our study, which are shown in Table 9.It can be seen that the elicited prior hyperparameters, in general, differ from the corresponding EB estimates (see Tables 6 and 9).Using these elicited parameters of Dirichlet prior, the final estimates of posterior cell probabilities, log odds ratios and interaction parameters were obtained for the dataset shown in Table 5.The estimated posterior means and the corresponding standard deviations of various cell probabilities are shown in Table 10.We also obtained the estimates of posterior cell probabilities for 2 × 2 × 2 setup but we are not showing those results, instead we shall focus on the quantities which are of more interest to the practitioners.The corresponding estimates for log odds ratio and other interactive parameters are shown in Table 11.It can be seen from Table 11 that θ GE is non-zero, which clearly conveys the message that genes and environment are associated with each other and cannot be considered independent.log(OR E ) is negative for both OC use and parity conveying that for non-susceptible subjects both the environmental components reduce the risk of ovarian cancer.Positive values of log(OR G ) suggest that patients having BRCA1/2 mutation are having high risk of disease in the control population.The estimated values of log(ψ) are showing positive association among genetic susceptibility and environmental exposure but in general these values are less than those obtained in the control population (refer the values of log(OR G ) in Table 11).Thus it can be concluded that both OC use and parity reduce the risk of ovarian cancer to a larger extent as compared to the categories of no-OC use and no-parity.It can be further seen that estimated variability along with the estimated posterior means convey that negative values for log(ψ) are also quite probable for both OC use and parity suggesting a further possibility of decrease in the values of log(ψ).
Table 12 provides the estimates of log odds ratios and interaction parameters for 2 × 2 × 4 setup.θ GE k , k = 1, ..., 3, is showing association for both OC use and parity discarding the assumption of independence between genetic susceptibility and environmental exposures.log(OR E k ), k = 1, ..., 3, is coming out to be negative suggesting that for non-susceptible subjects both OC use and parity reduce the risk of ovarian cancer and this risk mostly decreases for higher levels of OC and/or parity.The estimated values corresponding to log(OR G ) are not shown in Table 12.These values were 2.8567(0.2883)and 2.8220(1.0933),respectively, for OC use and parity, where bracketed values correspond to the estimated standard deviations.Therefore, it can be concluded that the patients with BRCA1/2 mutation are having high risk of developing ovarian cancer in control population.The log(ψ k ), k = 1, 2, 3, at different levels of environmental exposure are shown in the form of posterior density estimates using boxplot representations (see Figure 1).Although the density estimates at all the three levels of environmental exposures are showing more or less similar behaviour, they nicely convey an important message on the development of the disease with increasing OC and parity.
It can be seen that as the OC use and parity increase, the risk of ovarian cancer decreases.One or two estimates appear exception probably because of the reasons explained earlier.4

Results for an extracted data of reduced size
The results were also obtained for an extracted data of size 50 given in Table 3.This was done to see the effect of reduced sample size on our analysis.It is important to mention that if we have very small sample sizes, we may come across a situation where most of the cell frequencies are zero especially for large values of l.We, therefore, did not consider a value of n smaller than 50.The estimated cell probabilities along with their standard deviations (in parentheses) are given in Table 13.It is to be noted here that the prior hyperparameters are same that are presented in Table 9 for 2 × 2 × 2 set up.Contrary to our expectation, we observed that the estimated cell probabilities were showing more or less similar trend that we observed earlier for 2 × 2 × 2 setup.
We finally calculated the posterior estimates of log odds ratio and other interactive parameters, which  are shown in Table 14.It can be seen that the results are certainly changed in terms of numerical values but the final messages are more or less same that we drew earlier based on 2 × 2 × 2 setup for large sample size.That is, θ GE is showing an association among genetic and environmental components for both OC use and parity indicating the need of going a step ahead for evaluating the multiplicative interaction parameter ψ.Similarly, for non-susceptible subjects both OC use and parity reduce the risk of ovarian cancer whereas the subjects unexposed to environmental factors but having BRCA1/2 mutation are having a higher risk of the disease.The estimated ψ is also showing a similar behaviour with values indicating that OC use and parity reduce the risk of ovarian cancer even in patients with BRCA1/2 mutation.This reduction is stronger for OC use, than for parity.

Conclusion
The paper provides a complete Bayesian analysis of a gene-environment problem where the main focus is on prior elicitation of hyperparameters.The approach based on subjective elicitation of prior hyperparameters and subsequently the use of such elicited prior in drawing the final posterior inferences can always be regarded in a true Bayesian spirit.It is, however, often the case that availability of expert is not always guaranteed and, as such, the appropriate specification of prior hyperparameters becomes difficult.The methodology given here recommends the use of past data or any small subset of the given data to achieve the process of elicitation.It is seen that the considered methodology is, in general, easy to implement and provides a systematic way to tackle an important issue of hyperparameter selection.To the best of our understanding, no such attempt was ever made on the specific application considered in the paper where it is perhaps always desired to refine the inferences.The future researches may consider defining such strategies in general where one has priors other than the conjugate Dirichlet distribution.Elicitation of the complete functional form of the prior distribution is another important direction where future researchers may proceed.

Figure 1 :
Figure 1: Posterior density estimates of interaction parameter log(ψ k ), k = 1, 2, 3 (from left to right) corresponding to OC use and parity.

Table 2 :
Classification of case-control data with respect to disease status, genetic susceptibility and environmental exposure

Table 3 :
Classification of extracted data (2 × 2 × 2) of size 50 with respect to disease status, genetic susceptibility and environmental exposure

Table 7 :
Quartile values and estimates of hyperparameters corresponding to re-parameterized Dirichlet model for 2 × 2 × 4 setup

Table 8 :
Quartile values and estimates of hyperparameters corresponding to re-parameterized Dirichlet model for 2 × 2 × 2 setup

Table 9 :
Estimates of original Dirichlet hyperparameters

Table 10 :
Estimated sample based posterior means and the corresponding standard deviations of different cell probabilities for 2 × 2 × 4 set up

Table 11 :
Posterior estimates of log odds ratio and other association parameters for 2 × 2 × 2 setup

Table 12 :
Posterior estimates of log odds ratio and other association parameters using elicited prior hyperparameters (2 × 2 × 4 setup)

Table 13 :
Estimated sample based posterior means of different cell probabilities for extracted data of size 50

Table 14 :
Posterior estimates of log odds ratio and other association parameters for extracted data of size 50