# A Survey Design for a Sensitive Binary Variable Correlated with Another Nonsensitive Binary Variable

Authors: Jun-Wu Yu; Yang Lu; Guo-Liang Tian
Journal: Journal of Probability and Statistics
 (2013)

---

## Abstract

Tian et al. (2007) introduced a so-called hidden sensitivity model for evaluating the association of two sensitive questions with binary outcomes. However, in practice, we sometimes need to assess the association between one sensitive binary variable (e.g., whether or not a drug user, the number of sex partner being⩽1 or >1, and so on) and one nonsensitive binary variable (e.g., good or poor health status, with or without cervical cancer, and so on). To address this issue, by sufficiently utilizing the information contained in the non-sensitive binary variable, in this paper, we propose a new survey scheme, called combination questionnaire design/model, which consists of a main questionnaire and a supplemental questionnaire. The introduction of the supplemental questionnaire which is indeed a design of direct questioning can effectively reduce the noncompliance behavior since more respondents will not be faced with the sensitive question. Likelihood-based inferences including maximum likelihood estimates via the expectation-maximization algorithm, asymptotic confidence intervals, and bootstrap confidence intervals of parameters of interest are derived. A likelihood ratio test is provided to test the association between the two binary random variables. Bayesian inferences are also discussed. Simulation studies are performed, and a cervical cancer data set in Atlanta is used to illustrate the proposed methods.


---

## Body

## 1. Introduction

Warner [1] introduced a randomized response technique to obtain truthful answers to questions with sensitive attributes. Using the Warner design, Kraemer [2] derived a bivariate correlation between an attribute with polytomous responses and an attribute with normally distributed responses. Fox and Tracy [3] derived estimation of the Pearson product moment correlation coefficient between two sensitive questions by assuming that randomized response observations can be treated as individual-level scores that are contaminated by random measurement error. Edgell et al. [4] considered the correlation between two sensitive questions using the unrelated question design or the additive constants design. Christofides [5] presented a randomized response technique with two randomization devices to estimate the proportion of individuals having two sensitive characteristics at the same time. Kim and Warde [6] considered a multinomial randomized response model which can handle untruthful responses. They also derived the Pearson product moment correlation estimator which may be used to quantify the linear relationship between two variables when multinomial response data are observed according to a randomized response procedure. However, all these randomized response procedures make use of one or two randomizing devices which (i) entail extra costs in both efficiency and complexity, (ii) increase the cognitive load of randomized response techniques, and (iii) allow for new sources of error, such as misunderstanding the randomized response procedures or cheating on the procedures [7].From the perspective of incomplete categorical data design, Tian et al. [8] proposed a nonrandomized response model (called the hidden sensitivity model) for assessing the association of two sensitive questions with binary outcomes. To protect respondents' privacy and to avoid the use of any randomization device, they utilized a non-sensitive question in the questionnaire to indirectly obtain respondents' answers to the two sensitive questions. In the hidden sensitivity model, they implicitly assumed that all respondents are willing to follow the design instructions. In other words, the noncompliance behavior will not occur. However, in practice, we sometimes need to assess the association between one sensitive binary variable (e.g., whether or not a drug user, the number of sex partner being ⩽1 or >1, and so on) and one nonsensitive binary variable (e.g., good or poor health status, with or without cervical cancer, and so on). To our knowledge, the survey design for addressing this issue and corresponding statistical analysis methods are not available. Although we could directly adopt the hidden sensitivity model, the information contained in the nonsensitive binary variable cannot be utilized in the design. Intuitively, such information can be used to enhance the degree of privacy protection, so that more respondents will not be faced with the sensitive question. The major objective of this paper is to propose a new survey design and to develop corresponding statistical methods for analyzing sensitive data collected by this technique.The rest of the paper is organized as follows. In Section2, without using any randomizing device, we propose a survey scheme, called combination questionnaire design/model, which consists of a main questionnaire and a supplemental questionnaire. Likelihood-based inferences including maximum likelihood estimates via an expectation–maximization (EM) algorithm, asymptotic confidence intervals, and bootstrap confidence intervals of parameters of interest are derived in Section 3. A likelihood ratio test is also provided to test the association between the two binary random variables. In Section 4, we discuss Bayesian inferences when prior information on parameters is available. In Section 5, two simulation studies are performed to compare the efficiency of the proposed combination questionnaire model with that of the existing hidden sensitivity model of Tianet al. [8] (i.e., the main questionnaire only). A cervical cancer data set in Atlanta is used in Section 6 to illustrate the proposed methods. A discussion and an appendix on the mode of a group Dirichlet density and a sampling method from it are also presented.


## 2. The Survey Design

Assume thatX is a sensitive binary random variable, Y is a non-sensitive binary random variable, and they are correlated. Let {X=1} denote the sensitive class (e.g., X=1 if a respondent is a drug user) and {X=0} denote the non-sensitive class (e.g., X=0 if a respondent is not a drug user). Furthermore, let both {Y=1} (e.g., Y=1 if a respondent receives at least some college training) and {Y=0} (e.g., Y=0 if a respondent graduates at most from some high school) be non-sensitive classes. Define θ=(θ1,…,θ4)⊤, where θ1=Pr(X=0,Y=0), θ2=Pr(X=1,Y=0), θ3=Pr(X=1,Y=1), and θ4=Pr(X=0,Y=1), then θ∈𝕋4, where
(1)𝕋n=^{(x1,…,xn)⊤:xi⩾0,i=1,…,n,∑i=1n‍xi=1}
denotes the n-dimensional closed simplex in ℝn. The objective is to make inferences on θ, θx=^Pr(X=1)=θ2+θ3, θy=^Pr(Y=1)=θ4+θ3 and the odds ratio ψ=θ1θ3/(θ2θ4).The survey scheme consists of a main questionnaire and a supplemental questionnaire. To design the main questionnaire which is to be assigned to group 1 withn respondents (n is specified by the investigators), we first introduce a non-sensitive question (say, QW) with three possible answers. Assume that W is a non-sensitive variate with trichotomous outcomes associated with the QW and W is independent of both X and Y. Define pi=Pr(W=i) for i=1,2,3. For example, let W=1(2,3) denote that a respondent was born in January–April (May–August, September–December), and thus we could assume that pi≈1/3. The main questionnaire is shown in Table 1, under which each respondent is asked to answer the non-sensitive question QW.Table 1
The main questionnaire for the combination questionnaire model.
Category
W
=
1
W
=
2
W
=
3
I:{X=0,Y=0}
Block 1: —
Block 2: —
Block 3: —
II:{X=1,Y=0}
Category  II: Please put a tick in Block 2
III:{X=1,Y=1}
Category  III: Please put a tick in Block 3
IV:{X=0,Y=1}
Block 4: —
Note: Only{X=1} is a sensitive class, while {X=0}, {Y=0}, and {Y=1} are nonsensitive classes.On the one hand, since Category I (i.e.,{X=0,Y=0}) and Category IV (i.e., {X=0,Y=1}) are non-sensitive to each respondent, it is reasonable to assume that a respondent is willing to provide his/her truthful answer by putting a tick in Block i (i=1,…,4) according to his/her true status. On the other hand, Category II (i.e., {X=1,Y=0}) and Category III (i.e., {X=1,Y=1}) are usually sensitive to respondents. In this case, if a respondent belongs to Category II (III), he/she is designed to put a tick in Block 2 (3).The supplemental questionnaire is designed as shown in Table2, under which m respondents (m is also specified by the investigators) in group 2 are asked to put a tick in Block 5 or Block 6 depending on their true status; that is, {Y=0} or {Y=1}. Since both the {Y=0} and {Y=1} are non-sensitive classes, the supplemental questionnaire is in fact a design of direct questioning. Therefore, we call this design the combination questionnaire model.Table 2
The supplemental questionnaire for the combination questionnaire model.
Category
{
Y
=
0
}
Please put a tick in Block 5: — if you belong to{Y=0}
{
Y
=
1
}
Please put a tick in Block 6: — if you belong to{Y=1}
Note: Both the{Y=0} and {Y=1} are nonsensitive classes.Table3 shows the cell probabilities {θi}i=14, the observed frequencies {ni}i=14 and the unobservable frequencies {Zi}i=13, for the main questionnaire. The observed frequency n2 is the sum of the frequency of respondents belonging to Block 2 and the frequency of those belonging to Category II. The observed frequency n3 is the sum of the frequency of respondents belonging to Block 3 and the frequency of those belonging to Category III. Note that n=∑i=14‍ni=∑i=13‍Zi+n4, we have Z1=n-(Z2+Z3)-n4. Thus, only Z2 and Z3 are unobservable. Table 4 shows the cell probabilities and observed counts for the supplemental questionnaire.Table 3
Cell probabilities, observed and unobservable counts for the main questionnaire.
Category
W
=
1
W
=
2
W
=
3
Total
I:{X=0,Y=0}
p
1
θ
1
p
2
θ
1
p
3
θ
1
θ
1(Z1)
II:{X=1,Y=0}
θ
2(Z2)
III:{X=1,Y=1}
θ
3(Z3)
IV:{X=0,Y=1}
θ
4(n4)
θ
4(n4)
Total
p
1
θ
1
(
n
1
)
p
2
θ
1
+
θ
2
(
n
2
)
p
3
θ
1
+
θ
3
(
n
3
)
1(n)
Note:n=∑i=14‍ni, Z1=n-(Z2+Z3)-n4, where (Z2,Z3) are unobservable.Table 4
Cell probabilities and observed counts for the supplemental questionnaire.
Category
Total
{
Y
=
0
}
Block 5: —
θ
1
+
θ
2
(
m
0
)
{
Y
=
1
}
Block 6: —
θ
3
+
θ
4
(
m
1
)
Total
1(m)
Note:m=m0+m1.Remark 1.
The design of the main questionnaire is similar to that of the hidden sensitivity model of Tian et al. [8], while the design of the supplemental questionnaire is indeed a design of direct questioning since both the {Y=0} and {Y=1} are non-sensitive classes. Table 1 shows that Categories II and III are two sensitive subclasses. Therefore, putting a tick in Block 2 or Block 3 implies that the respondent could be suspected with the sensitive attribute. Let θ1=⋯=θ4≈0.25 (see Table 3); then around half of the, say, 2n respondents will be suspected with the sensitive attribute if only the main questionnaire is employed. However, besides the main questionnaire (with n respondents), if the supplemental questionnaire (see Table 2) with m=n respondents is also used, then only half of the n respondents will be suspected with the sensitive attribute. In other words, in the proposed combination questionnaire model, the information of the non-sensitive binary variable Y can be used to enhance the degree of privacy protection, so that more respondents will not be faced with the sensitive question. This is why we introduce the supplemental questionnaire besides the main questionnaire.Remark 2.
In practice, to simplify the design itself, we suggest that both the sample sizen in the main questionnaire and the sample size m in the supplemental questionnaire should be fixed in advance rather than the total sample size N=n+m being fixed. In this way, survey data can be collected in two independent groups, resulting in a relatively simpler statistical analysis. In addition, the interviewees are randomly assigned to either the group 1 or group 2.


## 3. Likelihood-Based Inferences

In this section,maximum likelihood estimates (MLEs) of the θ and the odds ratio ψ are derived by using the EM algorithm. In addition, asymptotic confidence intervals and the bootstrap confidence intervals of an arbitrary function of θ are also provided. Finally, a likelihood ratio test is presented for testing the association between the two binary random variables.


### 3.1. MLEs via the EM Algorithm

A total ofN=n+m respondents are classified into two groups by a randomization approach such that n respondents answer the questions in the main questionnaire and m respondents answer the questions in the supplemental questionnaire. Let Yobs,M={n;n1,…,n4} denote the observed counts collected in the main questionnaire (see Table 3), where ∑i=14‍ni=n. The likelihood function of θ based on Yobs,M is
(2)L(θ∣Yobs,M)=(nn1,n2,n3,n4)×(p1θ1)n1{∏i=23(piθ1+θi)ni}θ4n4.
Let Yobs,S={m;m0,m1} denote the observed counts gathered in the supplemental questionnaire (see Table 4), where m0+m1=m. The likelihood function of θ based on Yobs,S is
(3)L(θ∣Yobs,S)=(mm0)(θ1+θ2)m0(θ3+θ4)m1.
Let Yobs={Yobs,M,Yobs,S}. Since Yobs,M and Yobs,S are independent, the observed-data likelihood function of θ∈𝕋4 is
(4)LCQ(θ∣Yobs,S)∝θ1n1{∏i=23(piθ1+θi)ni}θ4n4×(θ1+θ2)m0(θ3+θ4)m1,
where the subscript “CQ” denotes the “combination questionnaire” model.Since the corresponding cell probabilities to the observed countsn2 and n3 in the group 1 are in the form of summation (i.e., piθ1+θi, i=2,3), we cannot obtain the explicit expressions of the MLEs of θ from the score equations of (4). By treating the observed counts n2 and n3 as incomplete data, we use the EM algorithm [9] to find the MLE θ^ of θ. The counts {Zi}i=13 in Table 3 can be viewed as missing data. Briefly, Z1,Z2,Z3, and n4 represent the counts of the respondents belonging to Categories I, II, III, and IV, respectively. Thus, we denote the latent data by Ymis={z2,z3} and the complete data by Ycom={Yobs,Ymis}. Note that all {pi} are known. Consequently, the complete-data likelihood function for θ is
(5)LCQ(θ∣YCOM)∝(∏i=13θizi)θ4n4×(θ1+θ2)m0(θ3+θ4)m1,
where z1=n-(z2+z3)-n4.By treating{θi}i=14 as random variables, we note that the complete-data likelihood function (5) has the density form of a grouped Dirichlet distribution [10]. Ng et al. [11] derived the mode of a grouped Dirichlet density with explicit expressions (see the appendix). Hence, from (A.4) and (A.5), the complete-data MLEs for θ are given by
(6)θ1=1-θ2-θ3-θ4,θ2=z2N(1+m0n-z3-n4),θ3=z3N(1+m1z3+n4),θ4=n4N(1+m1z3+n4).
Given Yobs and θ, Zi follows the binomial distribution with parameters ni and θi/(piθ1+θi); that is,
(7)Zi∣(Yobs,θ)~Binomial(ni,θipiθ1+θi),i=2,3.
Therefore, the E-step of the EM algorithm computes the following conditional expectations:
(8)E(Zi∣Yobs,θ)=niθipiθ1+θi,i=2,3,
and the M-step updates (6) by replacing z2 and z3 with previous conditional expectations.Remark 3.
Based on the observed-data likelihood function (4), we could use the Newton-Raphson algorithm to find the MLEs of θ. However, it is well known that the Newton-Raphson algorithm does not necessarily increase the log likelihood, leading even to divergence sometimes [12, page 172]. In addition, the Newton–Raphson algorithm is sensitive to the initial values. One advantage of using the EM algorithm in the current situation is that both the E- and M-step have closed-form expressions. More importantly, the EM algorithm and the data augmentation algorithm of Tanner and Wong [13] share the same data augmentation structure in the Bayesian settings (see Section 4 for more details).


### 3.2. Asymptotic Confidence Intervals

Letθ-4=(θ1,θ2,θ3)⊤. The asymptotic variance-covariance matrix of the MLE θ^-4 is then given by Iobs-1(θ^-4), where
(9)Iobs(θ-4)=-∂2ℓCQ(θ∣Yobs)∂θ-4∂θ-4⊤
denotes the observed information matrix and ℓCQ(θ∣Yobs)=logLCQ(θ∣Yobs) is the observed-data log-likelihood function. From (4), we have
(10)ℓCQ(θ∣Yobs)=n1logθ1+n2log(p2θ1+θ2)+n3log(p3θ1+θ3)+n4log(1-θ1-θ2-θ3)+m0log(θ1+θ2)+m1log(1-θ1-θ2).
It is easy to show that
(11)∂ℓCQ(θ∣Yobs)∂θ1=n1θ1-n4θ4+∑i=23‍nipipiθ1+θi+m0θ1+θ2-m11-θ1-θ2,∂ℓCQ(θ∣Yobs)∂θ2=-n4θ4+n2p2θ1+θ2+m0θ1+θ2-m11-θ1-θ2,∂ℓCQ(θ∣Yobs)∂θ3=-n4θ4+n3p3θ1+θ3,-∂2ℓCQ(θ∣Yobs)∂θ12=n1θ12+n4θ42+∑i=23‍nipi2(piθ1+θi)2+ϕ-∂2ℓCQ(θ∣Yobs)∂θ22=n4θ42+n2(p2θ1+θ2)2+ϕ-∂2ℓCQ(θ∣Yobs)∂θ32=n4θ42+n3(p3θ1+θ3)2-∂2ℓCQ(θ∣Yobs)∂θ1∂θ2=n4θ42+n2p2(p2θ1+θ2)2+ϕ-∂2ℓCQ(θ∣Yobs)∂θ1∂θ3=n4θ42+n3p3(p3θ1+θ3)2-∂2ℓCQ(θ∣Yobs)∂θ2∂θ3=n4θ42,
where
(12)ϕ=m0(θ1+θ2)2+m1(1-θ1-θ2)2.
Hence, the observed information matrix can be expressed as
(13)Iobs(θ-4)=diag(n1θ12,0,0)+n4θ42×1313⊤+A,
where
(14)A=(∑i=23nipi2(piθ1+θi)2+ϕ,n2p2(p2θ1+θ2)2+ϕ,n3p3(p3θ1+θ3)2n2p2(p2θ1+θ2)2+ϕ,n2(p2θ1+θ2)2+ϕ,0n3p3(p3θ1+θ3)2,0,n3(p3θ1+θ3)2).Letse(θ^i) denote the standard error of θ^i for i=1,2,3. Note that se(θ^i) can be estimated by the square root of the ith diagonal element of Iobs-1(θ^-4). We denote the estimated value of se(θ^i) by se^(θ^i). Thus, a 95% normal-based asymptotic confidence interval for θi can be constructed as
(15)[θ^i-1.96×se^(θ^i),θ^i+1.96×se^(θ^i)],i=1,2,3.Letϑ=h(θ-4) be an arbitrary differentiable function of θ-4. For example, θ4=1-∑i=13‍θi and the odds ratio ψ=θ1θ3/(θ2θ4). The delta method (e.g., [14, page 34]) can be used to approximate the standard error of ϑ^=h(θ^-4), and a 95% normal-based asymptotic confidence interval for ϑ is given by
(16)[ϑ^-1.96×se^(ϑ^),ϑ^+1.96×se^(ϑ^)],
where
(17)se^(ϑ^)={(∂ϑ∂θ-4)⊤Iobs-1(θ-4)(∂ϑ∂θ-4)|θ-4=θ^-4}1/2.


### 3.3. Bootstrap Confidence Intervals

When the normal-based asymptotic confidence interval like (15) is beyond the low bound zero or the upper bound one, the bootstrap approach [15] can be used to construct the bootstrap confidence interval of ϑ=h(θ-4). Based on the obtained MLE θ^, we independently generate
(18)(n1*,…,n4*)⊤~Multinomial(n;p1θ^1,p2θ^1+θ^2,p3θ^1+θ^3,θ^4),(19)m0*~Binomial(m,θ^1+θ^2).
Having obtained Yobs,M*={n;n1*,…,n4*} and Yobs,S*={m;m0*,m1*}, where m1*=m-m0*, we can calculate the bootstrap replication ϑ^*=h(θ-4*) based on Yobs*={Yobs,M*,Yobs,S*} via the EM algorithm specified by (6) and (8). Independently repeating this process G times, we obtain G bootstrap replications {ϑ^g*}g=1G. Consequently, a (1-α)100% bootstrap confidence interval for ϑ is given by
(20)[ϑ^L,ϑ^U],
where ϑ^L and ϑ^U are the 100(α/2) and 100(1-α/2) percentiles of {ϑ^g*}g=1G, respectively.


### 3.4. The Likelihood Ratio Test for Testing Association

The likelihood ratio statistic can be used to test whether the two binary random variablesX and Y are independent/correlated. The corresponding null and alternative hypotheses are [16, page 45]
(21)H0:ψ=1againstH1:ψ≠1.
The likelihood ratio statistic is defined by
(22)Λ=-2{ℓCQ(θ^R∣Yobs)-ℓCQ(θ^∣Yobs)},
where θ^R denotes the restricted MLE of θ under H0, θ^ denotes the MLE of θ, which can be obtained by the EM algorithm specified by (6) and (8), and ℓCQ(θ∣Yobs)=logLCQ(θ∣Yobs).To find the restricted MLEθ^R, we also employ the EM algorithm. Under H0:θ1θ3=θ2θ4, we have
(23)θ1=(1-θx)(1-θy),θ2=θx(1-θy),θ3=θxθy,θ4=(1-θx)θy.
In other words, under H0 we only have two free parameters θx and θy. Having obtained the restricted MLEs θ^x,R and θ^y,R, we can compute the restricted MLE θ^R=(θ^1,R,…,θ^4,R)⊤ from (23) by
(24)θ^1,R=(1-θ^x,R)(1-θ^y,R),θ^2,R=θ^x,R(1-θ^y,R),θ^3,R=θ^x,Rθ^y,R,θ^4,R=(1-θ^x,R)θ^y,R.In what follows, we consider the computation of the restricted MLEsθ^x,R and θ^y,R. Now, the complete-data likelihood function (5) becomes
(25)LCQ(θx,θy∣Ycom,H0)∝{(1-θx)(1-θy)}z1{θx(1-θy)}z2×(θxθy)z3{(1-θx)θy}n4(1-θy)m0θym1=θxz2+z3(1-θx)z1+n4θyz3+n4+m1(1-θy)z1+z2+m0
so that the restricted MLEs of θx and θy based on the complete-data are given by
(26)θx=z2+z3n,θy=z3+n4+m1N,
respectively. Thus, the M-step of the EM algorithm calculates (26), and the E-step computes the conditional expectations given in (8), where {θi} are defined in (23). Finally, under H0, Λ asymptotically follows chi-squared distribution with one degree of freedom.


## 3.1. MLEs via the EM Algorithm

A total ofN=n+m respondents are classified into two groups by a randomization approach such that n respondents answer the questions in the main questionnaire and m respondents answer the questions in the supplemental questionnaire. Let Yobs,M={n;n1,…,n4} denote the observed counts collected in the main questionnaire (see Table 3), where ∑i=14‍ni=n. The likelihood function of θ based on Yobs,M is
(2)L(θ∣Yobs,M)=(nn1,n2,n3,n4)×(p1θ1)n1{∏i=23(piθ1+θi)ni}θ4n4.
Let Yobs,S={m;m0,m1} denote the observed counts gathered in the supplemental questionnaire (see Table 4), where m0+m1=m. The likelihood function of θ based on Yobs,S is
(3)L(θ∣Yobs,S)=(mm0)(θ1+θ2)m0(θ3+θ4)m1.
Let Yobs={Yobs,M,Yobs,S}. Since Yobs,M and Yobs,S are independent, the observed-data likelihood function of θ∈𝕋4 is
(4)LCQ(θ∣Yobs,S)∝θ1n1{∏i=23(piθ1+θi)ni}θ4n4×(θ1+θ2)m0(θ3+θ4)m1,
where the subscript “CQ” denotes the “combination questionnaire” model.Since the corresponding cell probabilities to the observed countsn2 and n3 in the group 1 are in the form of summation (i.e., piθ1+θi, i=2,3), we cannot obtain the explicit expressions of the MLEs of θ from the score equations of (4). By treating the observed counts n2 and n3 as incomplete data, we use the EM algorithm [9] to find the MLE θ^ of θ. The counts {Zi}i=13 in Table 3 can be viewed as missing data. Briefly, Z1,Z2,Z3, and n4 represent the counts of the respondents belonging to Categories I, II, III, and IV, respectively. Thus, we denote the latent data by Ymis={z2,z3} and the complete data by Ycom={Yobs,Ymis}. Note that all {pi} are known. Consequently, the complete-data likelihood function for θ is
(5)LCQ(θ∣YCOM)∝(∏i=13θizi)θ4n4×(θ1+θ2)m0(θ3+θ4)m1,
where z1=n-(z2+z3)-n4.By treating{θi}i=14 as random variables, we note that the complete-data likelihood function (5) has the density form of a grouped Dirichlet distribution [10]. Ng et al. [11] derived the mode of a grouped Dirichlet density with explicit expressions (see the appendix). Hence, from (A.4) and (A.5), the complete-data MLEs for θ are given by
(6)θ1=1-θ2-θ3-θ4,θ2=z2N(1+m0n-z3-n4),θ3=z3N(1+m1z3+n4),θ4=n4N(1+m1z3+n4).
Given Yobs and θ, Zi follows the binomial distribution with parameters ni and θi/(piθ1+θi); that is,
(7)Zi∣(Yobs,θ)~Binomial(ni,θipiθ1+θi),i=2,3.
Therefore, the E-step of the EM algorithm computes the following conditional expectations:
(8)E(Zi∣Yobs,θ)=niθipiθ1+θi,i=2,3,
and the M-step updates (6) by replacing z2 and z3 with previous conditional expectations.Remark 3.
Based on the observed-data likelihood function (4), we could use the Newton-Raphson algorithm to find the MLEs of θ. However, it is well known that the Newton-Raphson algorithm does not necessarily increase the log likelihood, leading even to divergence sometimes [12, page 172]. In addition, the Newton–Raphson algorithm is sensitive to the initial values. One advantage of using the EM algorithm in the current situation is that both the E- and M-step have closed-form expressions. More importantly, the EM algorithm and the data augmentation algorithm of Tanner and Wong [13] share the same data augmentation structure in the Bayesian settings (see Section 4 for more details).


## 3.2. Asymptotic Confidence Intervals

Letθ-4=(θ1,θ2,θ3)⊤. The asymptotic variance-covariance matrix of the MLE θ^-4 is then given by Iobs-1(θ^-4), where
(9)Iobs(θ-4)=-∂2ℓCQ(θ∣Yobs)∂θ-4∂θ-4⊤
denotes the observed information matrix and ℓCQ(θ∣Yobs)=logLCQ(θ∣Yobs) is the observed-data log-likelihood function. From (4), we have
(10)ℓCQ(θ∣Yobs)=n1logθ1+n2log(p2θ1+θ2)+n3log(p3θ1+θ3)+n4log(1-θ1-θ2-θ3)+m0log(θ1+θ2)+m1log(1-θ1-θ2).
It is easy to show that
(11)∂ℓCQ(θ∣Yobs)∂θ1=n1θ1-n4θ4+∑i=23‍nipipiθ1+θi+m0θ1+θ2-m11-θ1-θ2,∂ℓCQ(θ∣Yobs)∂θ2=-n4θ4+n2p2θ1+θ2+m0θ1+θ2-m11-θ1-θ2,∂ℓCQ(θ∣Yobs)∂θ3=-n4θ4+n3p3θ1+θ3,-∂2ℓCQ(θ∣Yobs)∂θ12=n1θ12+n4θ42+∑i=23‍nipi2(piθ1+θi)2+ϕ-∂2ℓCQ(θ∣Yobs)∂θ22=n4θ42+n2(p2θ1+θ2)2+ϕ-∂2ℓCQ(θ∣Yobs)∂θ32=n4θ42+n3(p3θ1+θ3)2-∂2ℓCQ(θ∣Yobs)∂θ1∂θ2=n4θ42+n2p2(p2θ1+θ2)2+ϕ-∂2ℓCQ(θ∣Yobs)∂θ1∂θ3=n4θ42+n3p3(p3θ1+θ3)2-∂2ℓCQ(θ∣Yobs)∂θ2∂θ3=n4θ42,
where
(12)ϕ=m0(θ1+θ2)2+m1(1-θ1-θ2)2.
Hence, the observed information matrix can be expressed as
(13)Iobs(θ-4)=diag(n1θ12,0,0)+n4θ42×1313⊤+A,
where
(14)A=(∑i=23nipi2(piθ1+θi)2+ϕ,n2p2(p2θ1+θ2)2+ϕ,n3p3(p3θ1+θ3)2n2p2(p2θ1+θ2)2+ϕ,n2(p2θ1+θ2)2+ϕ,0n3p3(p3θ1+θ3)2,0,n3(p3θ1+θ3)2).Letse(θ^i) denote the standard error of θ^i for i=1,2,3. Note that se(θ^i) can be estimated by the square root of the ith diagonal element of Iobs-1(θ^-4). We denote the estimated value of se(θ^i) by se^(θ^i). Thus, a 95% normal-based asymptotic confidence interval for θi can be constructed as
(15)[θ^i-1.96×se^(θ^i),θ^i+1.96×se^(θ^i)],i=1,2,3.Letϑ=h(θ-4) be an arbitrary differentiable function of θ-4. For example, θ4=1-∑i=13‍θi and the odds ratio ψ=θ1θ3/(θ2θ4). The delta method (e.g., [14, page 34]) can be used to approximate the standard error of ϑ^=h(θ^-4), and a 95% normal-based asymptotic confidence interval for ϑ is given by
(16)[ϑ^-1.96×se^(ϑ^),ϑ^+1.96×se^(ϑ^)],
where
(17)se^(ϑ^)={(∂ϑ∂θ-4)⊤Iobs-1(θ-4)(∂ϑ∂θ-4)|θ-4=θ^-4}1/2.


## 3.3. Bootstrap Confidence Intervals

When the normal-based asymptotic confidence interval like (15) is beyond the low bound zero or the upper bound one, the bootstrap approach [15] can be used to construct the bootstrap confidence interval of ϑ=h(θ-4). Based on the obtained MLE θ^, we independently generate
(18)(n1*,…,n4*)⊤~Multinomial(n;p1θ^1,p2θ^1+θ^2,p3θ^1+θ^3,θ^4),(19)m0*~Binomial(m,θ^1+θ^2).
Having obtained Yobs,M*={n;n1*,…,n4*} and Yobs,S*={m;m0*,m1*}, where m1*=m-m0*, we can calculate the bootstrap replication ϑ^*=h(θ-4*) based on Yobs*={Yobs,M*,Yobs,S*} via the EM algorithm specified by (6) and (8). Independently repeating this process G times, we obtain G bootstrap replications {ϑ^g*}g=1G. Consequently, a (1-α)100% bootstrap confidence interval for ϑ is given by
(20)[ϑ^L,ϑ^U],
where ϑ^L and ϑ^U are the 100(α/2) and 100(1-α/2) percentiles of {ϑ^g*}g=1G, respectively.


## 3.4. The Likelihood Ratio Test for Testing Association

The likelihood ratio statistic can be used to test whether the two binary random variablesX and Y are independent/correlated. The corresponding null and alternative hypotheses are [16, page 45]
(21)H0:ψ=1againstH1:ψ≠1.
The likelihood ratio statistic is defined by
(22)Λ=-2{ℓCQ(θ^R∣Yobs)-ℓCQ(θ^∣Yobs)},
where θ^R denotes the restricted MLE of θ under H0, θ^ denotes the MLE of θ, which can be obtained by the EM algorithm specified by (6) and (8), and ℓCQ(θ∣Yobs)=logLCQ(θ∣Yobs).To find the restricted MLEθ^R, we also employ the EM algorithm. Under H0:θ1θ3=θ2θ4, we have
(23)θ1=(1-θx)(1-θy),θ2=θx(1-θy),θ3=θxθy,θ4=(1-θx)θy.
In other words, under H0 we only have two free parameters θx and θy. Having obtained the restricted MLEs θ^x,R and θ^y,R, we can compute the restricted MLE θ^R=(θ^1,R,…,θ^4,R)⊤ from (23) by
(24)θ^1,R=(1-θ^x,R)(1-θ^y,R),θ^2,R=θ^x,R(1-θ^y,R),θ^3,R=θ^x,Rθ^y,R,θ^4,R=(1-θ^x,R)θ^y,R.In what follows, we consider the computation of the restricted MLEsθ^x,R and θ^y,R. Now, the complete-data likelihood function (5) becomes
(25)LCQ(θx,θy∣Ycom,H0)∝{(1-θx)(1-θy)}z1{θx(1-θy)}z2×(θxθy)z3{(1-θx)θy}n4(1-θy)m0θym1=θxz2+z3(1-θx)z1+n4θyz3+n4+m1(1-θy)z1+z2+m0
so that the restricted MLEs of θx and θy based on the complete-data are given by
(26)θx=z2+z3n,θy=z3+n4+m1N,
respectively. Thus, the M-step of the EM algorithm calculates (26), and the E-step computes the conditional expectations given in (8), where {θi} are defined in (23). Finally, under H0, Λ asymptotically follows chi-squared distribution with one degree of freedom.


## 4. Bayesian Inferences

To derive the posterior mode ofθ, we employ the EM algorithm again. The latent data Ymis={z2,z3} are the same as those in Section 3.1. Based on the complete-data likelihood function (5), if the Dirichlet distribution Dirichlet (a1,…,a4) is adopted as the prior distribution of θ, then the complete-data posterior distribution is a grouped Dirichlet (GD) distribution with the formal definition given by (A.2); that is,
(27)f(θ∣Yobs,Ymis)=GD4,2,2(θ-4∣a,b),
where a=(a1+z1,a2+z2,a3+z3,a4+n4)⊤, b=(m0,m1)⊤, and z1=n-(z2+z3)-n4. The conditional predictive distribution is
(28)f(Ymis∣Yobs,θ)=∏i=23Binomial(zi|ni,θipiθ1+θi).Therefore, theM-step of the EM algorithm is to calculate the complete-data posterior mode:
(29)θ1=1-θ2-θ3-θ4,θ2=a2-1+z2a+-4+N(1+m0a1+a2-2+n-z3-n4),θ3=a3-1+z3a+-4+N(1+m1a3+a4-2+z3+n4),θ4=a4-1+n4a+-4+N(1+m1a3+a4-2+z3+n4),
where a+=∑i=14‍ai, and the E-step is to replace {zi}i=23 by the conditional expectations given by (8).In addition, based on (27) and (28), the data augmentation algorithm of Tanner and Wong [13] can be used to generate posterior samples of θ. A sampling method from (27) is given in the appendix.


## 5. Simulation Studies

In this section, two simulation studies are conducted to compare the efficiency of the proposed combination questionnaire model with that of the hidden sensitivity model of Tian et al. [8] (i.e., the main questionnaire model only), where (p1,p2,p3) are assumed to be (1/3,1/3,1/3) and pi=^Pr(W=i) for i=1,2,3. In the first simulated example, let the total sample size in the combination questionnaire model be the same as the sample size in the hidden sensitivity model. In the second example, the sample size for the main questionnaire in the combination questionnaire model is assumed to be equal to the sample size in the hidden sensitivity model.In the first simulated example, let a total ofN=n+m=50+50=100 participants be interviewed by using the combination questionnaire model. The true values of {θi}i=14 are listed in the second column of Table 5. We first generate
(30)(n1,…,n4)⊤~Multinomial(n;p1θ1,p2θ1+θ2,p3θ1+θ3,θ4),
and m0~Binomial(m,θ1+θ2) so that m1=m-m0. The EM algorithm (6) and (8) is used to calculate the MLEs of {θi}i=14. We repeated this experiment 1000 times. The average MLEs of {θi}i=14 are reported in the third column of Table 5. The corresponding bias, variance, and mean square error (MSE) are displayed in the fourth, fifth, and sixth columns of Table 5. Next, let N=n+m=100+0=100 participants be interviewed by using the hidden sensitivity model (i.e., the main questionnaire only). The corresponding results are reported in the last four columns of Table 5.Table 5
Comparison of the efficiency of the proposed combination questionnaire model with that of the hidden sensitivity model of Tian  et al.  [8] (i.e., the main questionnaire model only) with (p1,p2,p3)=(1/3,1/3,1/3).
Parameter
True value
The combination questionnaire model
The hidden sensitivity model
(N=n+m=50+50)
(N=n+m=100+0)
MLE
Bias
Variance
MSE
MLE
Bias
Variance
MSE
θ
1
0.10
0.1013
0.0013
0.0113
0.0113
0.1010
0.0010
0.0064
0.0064
θ
2
0.50
0.4985
-
0.0015
0.0089
0.0089
0.4975
-
0.0025
0.0041
0.0041
θ
3
0.20
0.2008
0.0008
0.0037
0.0037
0.2001
0.0001
0.0029
0.0029
θ
4
0.20
0.1994
-
0.0006
0.0026
0.0026
0.2014
0.0014
0.0016
0.0016
θ
1
0.20
0.1957
-
0.0043
0.0090
0.0090
0.1985
-
0.0015
0.0055
0.0055
θ
2
0.30
0.2999
-
0.0001
0.0070
0.0070
0.2999
-
0.0001
0.0034
0.0034
θ
3
0.40
0.4038
0.0038
0.0040
0.0040
0.4025
0.0025
0.0037
0.0037
θ
4
0.10
0.1006
0.0006
0.0017
0.0017
0.0991
-
0.0009
0.0009
0.0009
θ
1
0.30
0.3021
0.0021
0.0107
0.0107
0.3026
0.0026
0.0081
0.0081
θ
2
0.25
0.2483
-
0.0017
0.0077
0.0077
0.2499
-
0.0001
0.0039
0.0039
θ
3
0.15
0.1474
-
0.0026
0.0042
0.0042
0.1476
-
0.0024
0.0032
0.0032
θ
4
0.30
0.3021
0.0021
0.0033
0.0033
0.2999
-
0.0001
0.0021
0.0021
Note:Bias(θ^i)=E(θ^i)-θi and MSE(θ^i)=[Bias(θ^i)]2+Var(θ^i).From Table5, we can see that both the MLEs of {θi}i=14 in the combination questionnaire model and the hidden sensitivity model are very close to their true values, while the MSEs of {θ^i}i=14 in the combination questionnaire model are slightly larger than those in the hidden sensitivity model. These numerical results are not surprising since in the hidden sensitivity model, Tian et al. [8] implicitly assumed that all respondents must strictly follow the design instructions. In other words, the noncompliance behavior will not occur. The introduction of the supplemental questionnaire in the combination questionnaire model can effectively reduce the non-compliance behavior since more respondents will not be faced with the sensitive question, while the cost for introducing such a supplemental questionnaire is that we definitely lose a little of efficiency.In the second simulated example, we assume that a total ofN=n+m=100+100=200 participants are interviewed by using the combination questionnaire model, while only N=n+m=100+0=100 are interviewed by using the hidden sensitivity model. We repeat this experiment 1000 times. The corresponding results are reported in Table 6.Table 6
Comparison of the efficiency of the proposed combination questionnaire model with that of the hidden sensitivity model of Tian   et al.  [8] (i.e., the main questionnaire model only) with (p1,p2,p3)=(1/3,1/3,1/3).
Parameter
True value
The combination questionnaire model
The hidden sensitivity model
(N=n+m=100+100)
(N=n+m=100+0)
MLE
Bias
Variance
MSE
MLE
Bias
Variance
MSE
θ
1
0.10
0.1003
0.0003
0.0030
0.0030
0.1013
0.0013
0.0064
0.0064
θ
2
0.50
0.5015
0.0015
0.0029
0.0029
0.5019
0.0019
0.0041
0.0041
θ
3
0.20
0.1986
-
0.0014
0.0017
0.0017
0.1969
-
0.0031
0.0028
0.0028
θ
4
0.20
0.1995
-
0.0005
0.0013
0.0013
0.1999
-
0.0001
0.0016
0.0016
θ
1
0.20
0.2007
0.0007
0.0041
0.0041
0.2049
0.0049
0.0057
0.0057
θ
2
0.30
0.2989
-
0.0011
0.0033
0.0033
0.2965
-
0.0035
0.0034
0.0034
θ
3
0.40
0.4007
0.0007
0.0020
0.0020
0.3992
-
0.0008
0.0037
0.0037
θ
4
0.10
0.0997
-
0.0003
0.0009
0.0009
0.0994
-
0.0006
0.0009
0.0009
θ
1
0.30
0.3032
0.0032
0.0054
0.0054
0.2965
-
0.0036
0.0079
0.0079
θ
2
0.25
0.2488
-
0.0012
0.0039
0.0039
0.2501
0.0001
0.0038
0.0038
θ
3
0.15
0.1481
-
0.0019
0.0021
0.0021
0.1531
0.0031
0.0032
0.0032
θ
4
0.30
0.3000
0.0000
0.0017
0.0017
0.3003
0.0003
0.0021
0.0021
Note:Bias(θ^i)=E(θ^i)-θi and MSE(θ^i)=[Bias(θ^i)]2+Var(θ^i).From Table6, we can see that the MSEs of {θ^i}i=14 in the combination questionnaire model are smaller than those in the hidden sensitivity model. In addition, by comparing the fifth columns in Tables 5 and 6, we can see that the precisions of {θ^i}i=14 for the proposed combination questionnaire model in the second simulated example are significantly improved when compared with those in the first simulated example.


## 6. Analyzing Cervical Cancer Data in Atlanta

Williamson and Haber [17] reported a study which examined the relationship between disease status of cervical cancer and the number of sex partners and other risk factors. Cases were 20–79-year-old women of Fulton or DeKalb county in Atlanta, Georgia. They were diagnosed and were ascertained to have invasive cervical cancer. Controls were randomly chosen from the same counties and the same age ranges. Table 7 gives the cross-classification of number of sex partners (“few, 0–3” or “many, ⩾4”, denoted by X=0 or X=1) and disease status (control or case, denoted by Y=0 or Y=1). Generally, a sizable proportion (13.5% in this example) of the responses would be missing because of the sensitive question about the number of sex partners in a telephone interview. The objective is to examine if association exists between the number of sex partners and disease status of cervical cancer.Table 7
Cervical cancer data from Williamson and Haber [17].
Number of sex partners
Disease status of cervical cancer
Y
=
0 (control)
Y
=
1 (case)
X
=
0  (few, 0–3)
165
(
r
1
,
θ
1
)
103
(
r
4
,
θ
4
)
X
=
1 (many, ⩾4)
164
(
r
2
,
θ
2
)
221
(
r
3
,
θ
3
)
Missing
43
(
r
12
,
θ
1
+
θ
2
)
59
(
r
34
,
θ
3
+
θ
4
)
Note: The observed counts and the corresponding cell probabilities are in parentheses.X is a sensitive binary variate and Y is a non-sensitive binary variate.For the purpose of illustration, we presume that{X=0} is non-sensitive although the number 0–3 of sex partners is somewhat sensitive for some respondents. To illustrate the proposed design and approaches, we let W=1(2,3) if a respondent was born in January–April (May–August, September–December). It is then reasonable to assume that pk=Pr(W=k)≈1/3 for k=1,2,3, and W is independent of the sensitive binary variate X and the the non-sensitive binary variate Y. For the ideal situation (i.e., no sampling errors), the observed counts from the main questionnaire as shown in Tables 1 and 3 would be n1=r1/3=55, n2=55+r2=219, n3=55+r3=276, n4=r4=103; that is,
(31)Yobs,M={n;n1,…,n4}={653;55,219,276,103}.
On the other hand, we can view the missing data in Table 7 as the observed counts from the supplemental questionnaire as shown in Table 2. From Table 4, we have m0=r12=43 and m1=r34=59; that is, Yobs,S={m;m0,m1}={102;43,59}. Therefore, we obtain the observed data Yobs={Yobs,M,Yobs,S}.


### 6.1. Likelihood-Based Inferences

Usingθ(0)=14/4 as the initial values, the EM algorithm in (6) and (8) converged in 29 iterations. The resultant MLEs of θ and ψ are given in the second column of Table 8. From (13) and (14), the asymptotic variance-covariance matrix of the MLEs θ^-4 is
(32)Iobs-1(θ^-4)=(0.00086302-0.000442447-0.000384452-0.000442450.000511253-0.000010026-0.00038445-0.0000100260.000505241).
The estimated standard errors of θ^i (i=1,2,3) are square roots of the main diagonal elements of the above matrix. From (17), the estimated standard errors of θ^4=1-θ^1-θ^2-θ^3 and ψ^=θ^1θ^3/(θ^2θ^4) are given by
(33)se^(θ^4)={13⊤Iobs-1(θ^-4)13}1/2,se^(ψ^)={α⊤Iobs-1(θ^-4)α}1/2,
respectively, where
(34)α=(θ^3(1-θ^2-θ^3)θ^2θ^42,θ^1θ^3(θ^2-θ^4)(θ^2θ^4)2,θ^1(1-θ^1-θ^2)θ^2θ^42)⊤.Table 8
MLEs and95% confidence intervals of θ and ψ.
Parameter
MLE
std
95% asymptotic
95% bootstrap
CI
CI
θ
1
0.23793
0.02937
[0.18036, 0.29551]
[0.18006, 0.29879]
θ
2
0.24916
0.02261
[0.20484, 0.29348]
[0.20375, 0.29359]
θ
3
0.35196
0.02247
[0.30790, 0.39601]
[0.30731, 0.39666]
θ
4
0.16095
0.01434
[0.13284, 0.18905]
[0.13400, 0.18908]
ψ
2.08830
0.43971
[1.22646, 2.95011]
[1.35361, 3.18593]These estimated standard errors are listed in the third column of Table8. From (15) and (16), we can obtain the 95% asymptotic confidence intervals of θ and ψ, which are showed in the fourth column of Table 8.Based on (18) and (19), we generate G= 10,000 bootstrap samples. The corresponding 95% bootstrap confidence intervals of θ and ψ are displayed in the last column of Table 8.To test the null hypothesisH0:ψ=1 against H1:ψ≠1, we need to obtain the restricted MLE θ^R. Using
(35)(θx(0),θy(0))⊤=(0.5,0.5)⊤,
as the initial values, the EM algorithm in (26) converged to
(36)θ^x,R=0.64807,θ^y,R=0.52948,
in 19 iterations. From (24), the restricted MLEs of θ are obtained as
(37)θ^R=(0.16559,0.30493,0.34314,0.18634)⊤.
The log-likelihood ratio statistic Λ is equal to 12.469 and the P-value is 0.0004137. Since this P-value is far less than 0.05, the H0 is rejected at the 0.05 level of significance. Thus, we can conclude that there is an association between sex partners and cervical cancer status based on the current data. This conclusion is identical to that from the two 95% confidence intervals of the odds ratio as shown in Table 8, where both confidence intervals exclude the value 1.


### 6.2. Bayesian Inferences

When Dirichlet(1,1,1,1) (i.e., the uniform distribution on 𝕋4) is adopted as the prior distribution of θ, the posterior modes of θare equal to the corresponding MLEs. Using θ(0)=14/4 as the initial values, we employ the data augmentation algorithm to generate 1,000,000 posterior samples of θ and discard the first half of the samples. The Bayesian estimates of θ and ψ are given in Table 9. Since the lower bound of the Bayesian credible interval of the ψ is larger than 1, we believe that there is an association between the number of sexual partners and cervical cancer status.Table 9
Posterior modes and estimates of parameters for the cervical cancer data.
Parameter
Posterior
Bayesian
Bayesian
95% Bayesian
mode
mean
std
credible interval
θ
1
0.23793
0.24025
0.02934
[0.18574, 0.30050]
θ
2
0.24916
0.24789
0.02251
[0.20365, 0.29192]
θ
3
0.35196
0.35029
0.02246
[0.30621, 0.39432]
θ
4
0.16095
0.16155
0.01431
[0.13448, 0.19042]
ψ
2.08830
2.14538
0.45838
[1.39399, 3.18678]The posterior densities of the{θi}i=14 and ψ estimated by a kernel density smoother are plotted in Figures 1 and 2.The posterior densities of the{θi}i=14 estimated by a kernel density smoother based on the second 500,000 posterior samples of θgenerated by the data augmentation algorithm when the prior distribution is Dirichlet (1,1,1,1). (a) The posterior density of θ1; (b) The posterior density of θ2; (c) The posterior density of θ3; (d) The posterior density of θ4.
(a)
(b)
(c)
(d)Figure 2
The posterior density of the odds ratioψ estimated by a kernel density smoother based on the second 500,000 posterior samples of θ generated by the data augmentation algorithm when the prior distribution is Dirichlet (1,1,1,1).


## 6.1. Likelihood-Based Inferences

Usingθ(0)=14/4 as the initial values, the EM algorithm in (6) and (8) converged in 29 iterations. The resultant MLEs of θ and ψ are given in the second column of Table 8. From (13) and (14), the asymptotic variance-covariance matrix of the MLEs θ^-4 is
(32)Iobs-1(θ^-4)=(0.00086302-0.000442447-0.000384452-0.000442450.000511253-0.000010026-0.00038445-0.0000100260.000505241).
The estimated standard errors of θ^i (i=1,2,3) are square roots of the main diagonal elements of the above matrix. From (17), the estimated standard errors of θ^4=1-θ^1-θ^2-θ^3 and ψ^=θ^1θ^3/(θ^2θ^4) are given by
(33)se^(θ^4)={13⊤Iobs-1(θ^-4)13}1/2,se^(ψ^)={α⊤Iobs-1(θ^-4)α}1/2,
respectively, where
(34)α=(θ^3(1-θ^2-θ^3)θ^2θ^42,θ^1θ^3(θ^2-θ^4)(θ^2θ^4)2,θ^1(1-θ^1-θ^2)θ^2θ^42)⊤.Table 8
MLEs and95% confidence intervals of θ and ψ.
Parameter
MLE
std
95% asymptotic
95% bootstrap
CI
CI
θ
1
0.23793
0.02937
[0.18036, 0.29551]
[0.18006, 0.29879]
θ
2
0.24916
0.02261
[0.20484, 0.29348]
[0.20375, 0.29359]
θ
3
0.35196
0.02247
[0.30790, 0.39601]
[0.30731, 0.39666]
θ
4
0.16095
0.01434
[0.13284, 0.18905]
[0.13400, 0.18908]
ψ
2.08830
0.43971
[1.22646, 2.95011]
[1.35361, 3.18593]These estimated standard errors are listed in the third column of Table8. From (15) and (16), we can obtain the 95% asymptotic confidence intervals of θ and ψ, which are showed in the fourth column of Table 8.Based on (18) and (19), we generate G= 10,000 bootstrap samples. The corresponding 95% bootstrap confidence intervals of θ and ψ are displayed in the last column of Table 8.To test the null hypothesisH0:ψ=1 against H1:ψ≠1, we need to obtain the restricted MLE θ^R. Using
(35)(θx(0),θy(0))⊤=(0.5,0.5)⊤,
as the initial values, the EM algorithm in (26) converged to
(36)θ^x,R=0.64807,θ^y,R=0.52948,
in 19 iterations. From (24), the restricted MLEs of θ are obtained as
(37)θ^R=(0.16559,0.30493,0.34314,0.18634)⊤.
The log-likelihood ratio statistic Λ is equal to 12.469 and the P-value is 0.0004137. Since this P-value is far less than 0.05, the H0 is rejected at the 0.05 level of significance. Thus, we can conclude that there is an association between sex partners and cervical cancer status based on the current data. This conclusion is identical to that from the two 95% confidence intervals of the odds ratio as shown in Table 8, where both confidence intervals exclude the value 1.


## 6.2. Bayesian Inferences

When Dirichlet(1,1,1,1) (i.e., the uniform distribution on 𝕋4) is adopted as the prior distribution of θ, the posterior modes of θare equal to the corresponding MLEs. Using θ(0)=14/4 as the initial values, we employ the data augmentation algorithm to generate 1,000,000 posterior samples of θ and discard the first half of the samples. The Bayesian estimates of θ and ψ are given in Table 9. Since the lower bound of the Bayesian credible interval of the ψ is larger than 1, we believe that there is an association between the number of sexual partners and cervical cancer status.Table 9
Posterior modes and estimates of parameters for the cervical cancer data.
Parameter
Posterior
Bayesian
Bayesian
95% Bayesian
mode
mean
std
credible interval
θ
1
0.23793
0.24025
0.02934
[0.18574, 0.30050]
θ
2
0.24916
0.24789
0.02251
[0.20365, 0.29192]
θ
3
0.35196
0.35029
0.02246
[0.30621, 0.39432]
θ
4
0.16095
0.16155
0.01431
[0.13448, 0.19042]
ψ
2.08830
2.14538
0.45838
[1.39399, 3.18678]The posterior densities of the{θi}i=14 and ψ estimated by a kernel density smoother are plotted in Figures 1 and 2.The posterior densities of the{θi}i=14 estimated by a kernel density smoother based on the second 500,000 posterior samples of θgenerated by the data augmentation algorithm when the prior distribution is Dirichlet (1,1,1,1). (a) The posterior density of θ1; (b) The posterior density of θ2; (c) The posterior density of θ3; (d) The posterior density of θ4.
(a)
(b)
(c)
(d)Figure 2
The posterior density of the odds ratioψ estimated by a kernel density smoother based on the second 500,000 posterior samples of θ generated by the data augmentation algorithm when the prior distribution is Dirichlet (1,1,1,1).


## 7. Discussion

In this paper, we develop a general framework of design and analysis for the combination questionnaire model, which consists of a main questionnaire and a supplemental questionnaire. In fact, the main questionnaire (see Table1) is a generalization of the nonrandomized triangular model [18] and is then a special case of the multicategory triangular model [19]. The supplemental questionnaire (see Table 2) is a design of direct questioning. The introduction of the supplemental questionnaire can effectively reduce the noncompliance behavior since more respondents will not be faced with the sensitive question, while the cost for introducing such a supplemental questionnaire is that we definitely lose a little of efficiency. The combination questionnaire model can be used to gather information to evaluate the association between one sensitive binary variable and one non-sensitive binary variable.We note that the proposed combination questionnaire model has one limitation in applications; that is, it cannot be applied to situation where two categories{X=0} and {X=1} are sensitive like income. For example, let X=0 if his/her annual income is $25,000 or less and X=1 if his/her annual income is more than $25,000. For such cases, it is worthwhile to develop new designs to address this issue. One way is to replace the main questionnaire in Table 1 by a four-category parallel model (see [20, Section 4.1]). The other way is to employ the parallel model [21] to collect the information on X and to employ the design of direct questioning to collect information on Y then we could use the logistic regression to estimate the odds ratio, which is one of our further researches.



---

*Source: 827048-2013-06-23.xml*