Prediction of a Sensitive Feature under Indirect Questioning via Warner’s Randomized Response Technique and Latent Class Model

Abstract

We investigate the association of a sensitive characteristic or latent variable with observed binary random variables by the randomized response (RR) technique of Warner in his publication (Warner, S.L. J. Am. Stat. Assoc. 1965, 60, 63–69) and a latent class model. First, an expectation-maximization (EM) algorithm is provided to easily estimate the parameters of the null and alternative/full models for the association between a sensitive characteristic and an observed categorical random variable under the RR design of Warner’s paper above. The likelihood ratio test (LRT) is utilized to identify observed categorical random variables that are significantly related to the sensitive trait. Another EM algorithm is then presented to estimate the parameters of a latent class model constructed through the sensitive attribute and the observed binary random variables that are obtained from dichotomizing observed categorical random variables selected from the above LRT. Finally, two classification criteria are conducted to predict an individual in the sensitive or non-sensitive group. The practicality of the proposed methodology is illustrated with an actual data set from a survey study of the sexuality of first-year students, except international students, at Feng Chia University in Taiwan in 2016.

1. Introduction

Questionnaire surveys have been widely used to gather data for studies in various research fields, including behavioral science, education, sociology, economics, psychology, bio-medicine, etc., via Google form, email, phone call, or face-to-face interview. The information gathered from these surveys is utilized to estimate, compare, and forecast unknown population proportions of sensitive characteristics of interest. However, it will be likely to misreport the proportion of individuals who self-report some sensitive characteristic, such as political opinions, domestic violence, discrimination, abortion, drug use, gender identity, sexual behavior, anti-social behavior, exam fraud, plagiarism, illegal income, tax evasion, and gambling, or to refuse to answer when asked directly a sensitive question. This can lead to errors in analysis results and mistakes in statistical inferences. For example, in a study conducted by [1] that involved direct interviews with unemployment benefit recipients, 75% of the survey participants who had engaged in welfare or unemployment benefits deception denied doing so. In addition, according to [2], if respondents were forced to answer directly sensitive questions, the percentage of non-heterosexual participants in a community is typically underestimated. Therefore, various indirect questioning techniques, such as randomized response (RR), unmatched count, and crosswise model techniques, have been proposed to lessen potential bias due to non-response and social desirability response and, as a result, improve the reliability of data obtained from responses to sensitive topics to estimate better the population proportion of people bearing a sensitive attribute. See, e.g., [3,4,5,6] for more details.

The RR technique (RRT), initially introduced by Warner [3], is one of the most famous indirect questioning techniques and has been widely used to collect sensitive data. The main idea of this technique is to safeguard respondents’ privacy by concealing their answers or without disclosing their actual status by using a random device, such as spinners, a deck of cards, dice, or random number generators. Specifically, Warner [3] suggested the related-question RRT via a sensitive question of interest and its complement. Respondents can feel more comfortable choosing which question to honestly answer “yes” or “no” based on the results of running a random device because the interviewer does not know precisely their responses. However, Warner’s RRT still has limitations, such as not working when the probability that the sensitive question is selected to answer is 0.5. Therefore, various extensions of Warner’s RRT have been suggested to overcome its limitations and enhance computing effectiveness. For example, Horvitz et al. [7] and Greenberg et al. [8] proposed an unrelated-question RR design in which the first one is a sensitive question, as in Warner’s design, and the second one is innocuous and independent of the sensitive question. Mangat and Singh [9] proposed a two-stage RR design that used two random devices in the procedure. Christofides [10] suggested a generalized RRT. Huang [11] applied the two-stage RRT to enhance the performance of Warner’s RRT. Tian and Tang [12] provided another classification for the RRTs. The related literature can also be found in [2,5,13,14,15,16,17,18,19,20,21].

Data of some auxiliary variables are also gathered under direct questions when collecting sensitive trait or behavior data. The effects of these auxiliary variables on the population proportion of a sensitive characteristic are also quite essential and mentioned in many studies. For instance, by using the RRT of [3], Maddala [22] estimated the effects of the auxiliary variables, e.g., sex, age, and place of residence, on drug use by utilizing a logistic regression model. Scheers and Dayton [23] obtained sensitive information through the RR designs of [3,8], respectively. They proposed using the logistic regression model to establish the influence model of the accompanying variables on the sensitive feature and to provide a method for estimating the regression parameters. Hsieh et al. [14,24] proposed semiparametric methods of estimating the parameters of a logistic regression model when values of some of the covariates on some subjects are missing at random to gather information about sensitive characteristics based on the RR design of [8]. Chang et al. [17] developed a covariate extension of the two-stage RR design of [11] by employing logistic regression to investigate the effects of two covariates of interest on the sensitive attribute and the truthful response to the directly asked sensitive question in the first stage for a sensitive trait respondent to estimate the probability of a sensitive feature respondent’s honest response in the first stage and the probability of a respondent with the sensitive feature based on the logistic regression parameter estimates. Furthermore, in some practical applications, in addition to using an RR design to collect data on sensitive characteristics or behaviors, the related scale data are also frequently gathered to use as auxiliary variables. For example, Lee et al. [19] studied students’ attitudes toward love, gender identity, and online dating experience through the survey study of the sexuality of first-year students, except international students, at Feng Chia University in Taiwan in 2016. They used the sum of scores of responses to six internet dating experience statements (statements 40–45) as the auxiliary variable for the sensitive issue, question 67: “Have you ever had a one-night stand through a dating site or mobile app?” in which response data were gathered by the generalized RRT of Christofides [10]. They utilized the proposed Bayesian estimation methods to estimate the parameters of the normal independent and dependent models for the association of the response to question 67 with the sum of scores of responses to the six internet dating experiences. Two Bayesian model selection criteria were proposed to choose one of the two models as an appropriate model to describe this association.

In another approach to sample survey research, several latent variable models, which are statistical models and relate a set of observed variables (so-called indicators/manifest variables) to a latent variable, have been widely used in various sciences, such as the social, economic, behavioral, and health sciences. The main types of latent variable models can be classified based on whether the manifest and latent variables are categorical or continuous as in (Table 7.1 [25], p. 178). See, e.g., [13,25,26] for more details. Some other extensions of latent variable models, where the manifest variables are treated as ordinal categorical variables, count variables, or other metrics, can be found in [27,28].

Because both the latent and indicators/manifest/observed variables are categorical, the latent class (LC) analysis (LCA) is recommended for applications. Lazarsfeld [29] first introduced the conceptual foundation of LCA, which was used as a tool for building typologies based on observed dichotomous variables, in research about the ethnocentrism of American soldiers during World War II [28,30]. Lazarsfeld and Henry [31] provided a thorough and in-depth conceptual and mathematical treatment of the LCA, but there was a lack of reliable parameter estimation methods. Since then, the LC model (LCM) has been improved and applied in various science fields. For example, Goodman [32,33] provided a relatively simple method, which was later shown to be closely related to the expectation-maximization (EM) algorithm [34], to obtain the maximum likelihood (ML) estimates of LCM parameters and goodness-of-fit methods for the model for fitting the observed data. Haberman [35] demonstrated the connection between LCMs and log-linear models for frequency tables with missing (unknown) cell counts. Dayton and Macready [36] proposed a new development by incorporating covariates into an LCM. Muthén and Shedden [37] improved the models that identified latent growth trajectory class membership in longitudinal data based on individual growth trajectories and were estimated by the EM algorithm. Stern et al. [38] applied LCA in studying the two main temperamental types of children: inhibited and uninhibited. LCMs have been used to investigate the initiation of substance use habits throughout adolescence, such as alcohol, caffeine, and tobacco [39]. Vermunt [40] provided an overview of applications of LCMs in social science research and an extension of the LCM to deal with nested data structures. Collins and Lanza [26] presented the methodology and applications of LCA for social, behavioral, and health sciences data. Nasiopoulou et al. [41] applied LCA to investigate the professional profiles of Swedish preschool teachers. LCA was used to determine the cause of occupational fatalities in the study of [42]. Wu et al. [43] used LCA to stratify the risk of incident diabetes in Chinese adults. For other studies about LCA and its applications, see, e.g., [26,44,45,46,47,48].

However, to the best of our knowledge, there has not been any research about applying the combination of the RRT and LCM to estimate the proportion of a sensitive or latent characteristic based on observed variables. Therefore, we are strongly motivated by the issue to introduce the proposed models, estimation methods, likelihood ratio test (LRT), and reality applications. This study first presents an EM algorithm for estimating the parameters of the null and alternative/full model for the association between the sensitive attribute and an observed categorical random variable in which the RRT of Warner [3] is used. Note that our approaches are different from those of Lee et al. [19] in which they used a Markov chain Monte Carlo estimation method and generalized RRT of [10]. The LRT is then applied to assess whether there is a difference in the distribution of the observed categorical random variable (i.e., auxiliary variable) between the sensitive and non-sensitive groups. Finally, a combination of the RRT of [3] and LCA is introduced in which the indicators/manifest/observed variables are binary. After collecting questions/statements significantly related to the sensitive characteristic from the results of the LRTs, an LCM is used to classify an individual in the sensitive or non-sensitive group.

Section 2.1 reviews the RR design of [3]. In Section 2.2, we introduce a model for the association between a sensitive variable and an observed categorical random variable, which is presented in [19], by using the RRT of [3]. The EM algorithm, called EM algorithm 1, is applied to estimate the parameters of this model, and the LRT is used to evaluate whether the observed categorical random variable is significantly associated with the sensitivity characteristic variable. Section 2.3 provides an LCM to relate observed dichotomous random variables, obtained from dichotomizing the significant observed categorical random variables, to a latent variable that is the sensitive attribute from the RR design of [3]. Another EM algorithm, called EM algorithm 2, is proposed to estimate the parameters of the LCM. A classification criterion is conducted to predict an individual in a sensitive or non-sensitive group in this section. In Section 3, data from the survey of the sexuality of freshmen, except international students, at Feng Chia University in Taiwan in 2016 are used to demonstrate practical applications of the proposed methodology. Another classification criterion is proposed to forecast students in the sensitive or non-sensitive group. Conclusions and some remarks are given in Section 4.

2. Models and Methods

2.1. Warner’s RR Design

The fundamental idea of the related-question RR design of [3] is to protect respondents’ privacy by concealing their responses via a random device, e.g., spinners and dice, a deck of cards, etc. For example, assume that each respondent receives randomly a deck of cards marked with question A or $\overline{A}$ as follows:

• $A:$ Have you ever had a one-night stand through a dating site or mobile app?
• $\overline{A}:$ Have you never had a one-night stand through a dating site or mobile app?

Respondents report a truthful “yes” or “no” response to the question they receive without showing the interviewer which question has been selected. The actual status of the respondents regarding whether or not they have ever had a one-night stand through a dating site or mobile app remains undisclosed, and, therefore, their privacy is protected because neither the interviewer nor the researcher is even aware of the question to which the released answer refers.

Let $\theta$ be the population proportion of persons bearing the sensitive trait having ever had a one-night stand through a dating site or mobile app, p the probability of cards marked with the question A, and $1-p$ the probability of cards marked with the question $\overline{A}$. Assume that p is known. The probability of answering “yes” is then

$$\begin{array}{cc}\hfill P(\mathrm{yes})& =p(A)P(\mathrm{yes}|A)+P(\overline{A})P(\mathrm{yes}|\overline{A})\hfill \\ \hfill & =\theta p+(1-\theta )(1-p)\hfill \\ \hfill & =\theta (2p-1)+1-p.\hfill \end{array}$$

Let $\widehat{\lambda }$ be the proportion of respondents who answer “yes”. Because $2p-1\ne 0$, the estimate of $\theta$ can be obtained as follows:

$$\begin{array}{c}\hfill {\widehat{\theta }}_w=\frac{\widehat{\lambda }-(1-p)}{2p-1}.\end{array}$$

However, using random devices may make a respondent report different answers when repeating the survey twice. Yu et al. [5] pointed out that this is a repetitive phenomenon and, hence, provided another design in which a respondent chooses question A or $\overline{A}$ to answer according to her/his characteristic. For example, if a respondent were born between August 11 and December 31, she/he chooses question $\overline{A}$ to answer and chooses question A to answer otherwise. Therefore, if the respondent repeated the survey twice, the results are the same. In addition, it can prevent respondents from rejecting interviews and providing false answers and, hence, improve estimation efficiency through this random birthday design to obtain data of responses to sensitive questions. In this study, for simplicity, we use the birthday design of [5] for the RR design of [3] by using the data in the study of [19] as reality analysis. By this approach, it is only necessary that specific explanations are written in the questionnaire, and students can then answer questions by themselves without using a random device. The details of specific explanations are given in Section 3.

2.2. Model for the Association between Sensitive Attribute and Categorical Variable

2.2.1. Model

Let $Y_i$, $i=1,2,\dots ,n$, denote whether the ith sample respondent has the sensitive attribute A in which $Y_i=1$ if yes and $Y_i=0$ otherwise. Assume that $Z_i$ denotes the ith respondent’s response to a direct inquiry, which is a categorical random variable with values $1,2,\dots ,B$ or a quantitative variable.

The model for the association between $Y_i$ and $Z_i$ is considered as follows:

$$\begin{array}{c}P(Y_i=y)={\theta }^y{(1-\theta )}^{1-y},\hfill \\ P(Z_i=z|Y_i=y)={\pi }_{yz},z=1,2,\dots ,B,y=0,1.\hfill \end{array}$$

(1)

From the above model, one can then obtain

$$\begin{array}{c}\hfill P(Z_i=z,Y_i=y)={\theta }^y{(1-\theta )}^{1-y}{\pi }_{yz},z=1,2,\dots ,B,y=0,1.\end{array}$$

(2)

Let $Y_i^0$ be the result of using this non-random answering design as given in Table 10 of [17]. Let $Y_i^0=1$ denote the answer of “A” to a sensitive question and $Y_i^0=0$ the answer of “B” to the sensitive question. Assume that p is the probability of respondents’ birthday between January 1 and August 10. $P(Z_i=z,Y_i^0=y^0)$ can then be expressed as follows:

$$\begin{array}{cc}\hfill P(Z_i=z,Y_i^0=1)& =P(Z_i=z,Y_i^0=1,Y_i=1)+P(Z_i=z,Y_i^0=1,Y_i=0)\hfill \\ \hfill & =P(Z_i=z,Y_i^0=1|Y_i=1)P(Y_i=1)+P(Z_i=z,Y_i^0=1|Y_i=0)P(Y_i=0)\hfill \\ \hfill & ={\pi }_{1z}p\theta +{\pi }_{0z}(1-p)(1-\theta ).\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill P(Z_i=z,Y_i^0=0)& =P(Z_i=z,Y_i^0=0,Y_i=1)+P(Z_i=z,Y_i^0=0,Y_i=0)\hfill \\ \hfill & =P(Z_i=z,Y_i^0=0|Y_i=1)P(Y_i=1)+P(Z_i=z,Y_i^0=0|Y_i=0)P(Y_i=0)\hfill \\ \hfill & ={\pi }_{1z}(1-p)\theta +{\pi }_{0z}p(1-\theta ).\hfill \end{array}$$

Let ${\mathit{y}}^0=(y_1^0,y_2^0,\dots ,y_n^0)$, $\mathit{z}=(z_1,z_2,\dots ,z_n)$, ${\mathbf{\pi }}_y=({\pi }_{y1},{\pi }_{y2},\dots ,{\pi }_{yB})$, $y=0,1$, and $\Theta =(\theta ,{\mathbf{\pi }}_0,{\mathbf{\pi }}_1)$. Because only $(Y_i^0,Z_i)$, $i=1,2,\dots ,n,$ can be observed, the observed-data likelihood function of $\Theta$ can then be written as

$$\begin{array}{cc}\hfill & L_{obs}(\Theta |{\mathit{y}}^0,\mathit{z})\hfill \\ \hfill =& \prod _{i=1}^n\left\{{\left[p\theta {\pi }_{1z_i}+(1-p)(1-\theta ){\pi }_{0z_i}\right]}^{I(y_i^0=1)}{\left[p(1-\theta ){\pi }_{0z_i}+(1-p)\theta {\pi }_{1z_i}\right]}^{I(y_i^0=0)}\right\}\hfill \\ \hfill =& \prod _{i=1}^n\prod _{j=1}^B{\left\{{\left[p\theta {\pi }_{1j}+(1-p)(1-\theta ){\pi }_{0j}\right]}^{I(y_i^0=1)}{\left[p(1-\theta ){\pi }_{0j}+(1-p)\theta {\pi }_{1j}\right]}^{I(y_i^0=0)}\right\}}^{I(z_i=j)}\hfill \\ \hfill =& \prod _{j=1}^B\left\{{\left[p\theta {\pi }_{1j}+(1-p)(1-\theta ){\pi }_{0j}\right]}^{{\sum }_{i=1}^{n}I(y_i^0=1,z_i=j)}{\left[p(1-\theta ){\pi }_{0j}+(1-p)\theta {\pi }_{1j}\right]}^{{\sum }_{i=1}^{n}I(y_i^0=0,z_i=j)}\right\}.\hfill \end{array}$$

Given an initial value for $\Theta$, the $\mathsf{R}$ function optim or nlminb can then be used to obtain the ML estimate of $\Theta$. However, in practice, when using the $\mathsf{R}$ function optim or nlminb to find the ML estimate of $\Theta$ for the RR design of [3], it may lead to encountering the problem of non-convergence, or the estimate is not in the parameter space of $\Theta$. Therefore, an EM algorithm is proposed to overcome this problem in the next section.

2.2.2. EM Algorithm 1

The EM algorithm developed by Dempster et al. [34] has been the most widely used iterative technique for computing ML estimates from incomplete data. The EM algorithm consists of two steps: the expectation (E)-step and maximization (M)-step. In the E-step, the latent or unobserved data are estimated by their expectation given the observed data and current parameter estimates. The M-step is used to maximize the expectation in the E-step to update estimates of unknown parameters.

Let $\mathit{y}=(y_1,y_2,\dots ,y_n)$. The complete-data likelihood function of $\Theta$ given $(\mathit{y},{\mathit{y}}^0,\mathit{z})$ is written as follows:

$$\begin{array}{cc}\hfill L_c(\Theta |\mathit{y},{\mathit{y}}^0,\mathit{z})& =\prod _{i=1}^n\left\{{\left[p\theta {\pi }_{1z_i}\right]}^{I(y_i^0=1,y_i=1)}{\left[(1-p)\theta {\pi }_{1z_i}\right]}^{I(y_i^0=0,y_i=1)}\right.\hfill \\ \hfill & \times \left.{\left[(1-p)(1-\theta ){\pi }_{0z_i}\right]}^{I(y_i^0=1,y_i=0)}{\left[p(1-\theta ){\pi }_{0z_i}\right]}^{I(y_i^0=0,y_i=0)}\right\}\hfill \\ \hfill & =\prod _{i=1}^n\prod _{j=1}^B\left\{{\left[p\theta {\pi }_{1j}\right]}^{y_i^0{y}_i}{\left[(1-p)\theta {\pi }_{1j}\right]}^{(1-y_i^0)y_i}\right.\hfill \\ \hfill & \times {\left.{\left[(1-p)(1-\theta ){\pi }_{0j}\right]}^{y_i^0(1-y_i)}{\left[p(1-\theta ){\pi }_{0j}\right]}^{(1-y_i^0)(1-y_i)}\right\}}^{I(z_i=j)}.\hfill \end{array}$$

The complete-data log-likelihood function of $\Theta$ given $(\mathit{y},{\mathit{y}}^0,\mathit{z})$ can then be expressed as follows:

$$\begin{array}{cc}\hfill {\ell }_c(\Theta |\mathit{y},{\mathit{y}}^0,\mathit{z})& =\sum _{j=1}^B\sum _{i=1}^n\left\{y_i^0{y}_{i}I(z_i=j)ln(p\theta {\pi }_{1j})+(1-y_i^0)y_{i}I(z_i=j)ln\left[(1-p)\theta {\pi }_{1j}\right]\right.\hfill \\ \hfill & +y_i^0(1-y_i)I(z_i=j)ln\left[(1-p)(1-\theta )\theta {\pi }_{0j}\right]\hfill \\ \hfill & \left.+(1-y_i^0)(1-y_i)I(z_i=j)ln\left[p(1-\theta ){\pi }_{0j}\right]\right\}.\hfill \end{array}$$

(3)

Note that although Y is a latent or unobserved variable, one can obtain $P(Y_i=y|Z_i=j,Y_i^0=y^0)$, $y=0,1$, $y^0=0,1$, $j=1,2,\dots ,B$, as follows:

$$\begin{array}{cc}\hfill P(Y_i=1|Z_i=j,Y_i^0=1)& ={\displaystyle \frac{p\theta {\pi }_{1j}}{p\theta {\pi }_{1j}+(1-p)(1-\theta ){\pi }_{0j}}},\hfill \\ \hfill P(Y_i=1|Z_i=j,Y_i^0=0)& ={\displaystyle \frac{(1-p)\theta {\pi }_{1j}}{(1-p)\theta {\pi }_{1j}+p(1-\theta ){\pi }_{0j}}},\hfill \\ \hfill P(Y_i=0|Z_i=j,Y_i^0=1)& ={\displaystyle \frac{(1-p)(1-\theta ){\pi }_{0j}}{p\theta {\pi }_{1j}+(1-p)(1-\theta ){\pi }_{0j}}},\hfill \\ \hfill P(Y_i=0|Z_i=j,Y_i^0=0)& ={\displaystyle \frac{p(1-\theta ){\pi }_{0j}}{(1-p)\theta {\pi }_{1j}+p(1-\theta ){\pi }_{0j}}}.\hfill \end{array}$$

Therefore, it can yield the following results:

$$\begin{array}{cc}\hfill E\left(Y_i|Z_i=j,Y_i^0=1\right)& =P(Y_i=1|Z_i=j,Y_i^0=1)\hfill \\ \hfill & ={\displaystyle \frac{p\theta {\pi }_{1j}}{p\theta {\pi }_{1j}+(1-p)(1-\theta ){\pi }_{0j}}},\hfill \end{array}$$

$$\begin{array}{cc}\hfill E\left(Y_i|Z_i=j,Y_i^0=0\right)& =P(Y_i=1|Z_i=j,Y_i^0=0)\hfill \\ \hfill & ={\displaystyle \frac{(1-p)\theta {\pi }_{1j}}{(1-p)\theta {\pi }_{1j}+p(1-\theta ){\pi }_{0j}}},\hfill \end{array}$$

$$\begin{array}{cc}\hfill E\left((1-Y_i)|Z_i=j,Y_i^0=1\right)& =P(Y_i=0|Z_i=j,Y_i^0=1)\hfill \\ \hfill & ={\displaystyle \frac{(1-p)(1-\theta ){\pi }_{0j}}{p\theta {\pi }_{1j}+(1-p)(1-\theta ){\pi }_{0j}}},\hfill \end{array}$$

$$\begin{array}{cc}\hfill E\left((1-Y_i)|Z_i=j,Y_i^0=0\right)& =P(Y_i=0|Z_i=j,Y_i^0=0)\hfill \\ \hfill & ={\displaystyle \frac{p(1-\theta ){\pi }_{0j}}{(1-p)\theta {\pi }_{1j}+p(1-\theta ){\pi }_{0j}}}.\hfill \end{array}$$

E-step: E-step is to take the expectation of ${\ell }_c(\Theta |\mathit{y},{\mathit{y}}^0,\mathit{z})$ in (3) with respect to the conditional distributions of the unobserved variables $Y_i$s given the current estimate of $\Theta$ and the observed data $({\mathit{y}}^0,\mathit{z})$. Let ${\widehat{\Theta }}^{(m)}=({\widehat{\theta }}^{(m)},{\widehat{\mathbf{\pi }}}_0^{(m)},{\widehat{\mathbf{\pi }}}_1^{(m)})$ be an estimate of $\Theta$ at the mth iteration, where ${\widehat{\mathbf{\pi }}}_y^{(m)}=({\widehat{\pi }}_{y1}^{(m)},{\widehat{\pi }}_{y2}^{(m)},\dots ,{\widehat{\pi }}_{yB}^{(m)})$, $y=0,1$. ${\widehat{\theta }}^{(0)}=max\{0.001,[{\overline{y}}^0-(1-p)]/(2p-1)\}$, where ${\overline{y}}^0=n^{-1}{\sum }_{i=1}^n{y}_i^0$, and ${\widehat{\pi }}_{1j}^{(0)}={\widehat{\pi }}_{0j}^{(0)}=n^{-1}{\sum }_{i=1}^{n}I(z_i=j)$ are initial values. Given the observed data $({\mathit{y}}^0,\mathit{z})$ and ${\widehat{\Theta }}^{(m)}$, by taking the expectation of ${\ell }_c(\Theta |\mathit{Y},{\mathit{Y}}^0,\mathit{Z})$ in (3), where ${\mathit{Y}}^0=(Y_1^0,Y_2^0,\dots ,Y_n^0)$, $\mathit{Y}=(Y_1,Y_2,\dots ,Y_n)$, and $\mathit{Z}=(Z_1,Z_2,\dots ,Z_n)$, the Q-function can be given as follows:

$$\begin{array}{cc}\hfill Q\left(\Theta |{\widehat{\Theta }}^{(m)}\right)& =E\left[{\ell }_c(\Theta |\mathit{Y},{\mathit{Y}}^0,\mathit{Z})|{\mathit{y}}^0,\mathit{z},{\widehat{\Theta }}^{(m)}\right]\hfill \\ \hfill & =\sum _{j=1}^B\left\{\left[\sum _{i=1}^n{y}_i^{0}I(z_i=j)E\left(Y_i|{\mathit{y}}^0,\mathit{z},{\widehat{\Theta }}^{(m)}\right)\right]ln(p\theta {\pi }_{1j})\right\}\hfill \\ \hfill & +\sum _{j=1}^B\left\{\left[\sum _{i=1}^n(1-y_i^0)I(z_i=j)E\left(Y_i|{\mathit{y}}^0,\mathit{z},{\widehat{\Theta }}^{(m)}\right)\right]ln\left[(1-p)\theta {\pi }_{1j}\right]\right\}\hfill \\ \hfill & +\sum _{j=1}^B\left\{\left[\sum _{i=1}^n{y}_i^{0}I(z_i=j)E\left(1-Y_i|{\mathit{y}}^0,\mathit{z},{\widehat{\Theta }}^{(m)}\right)\right]ln\left[(1-p)(1-\theta )\theta {\pi }_{0j}\right]\right\}\hfill \\ \hfill & +\sum _{j=1}^B\left\{\left[\sum _{i=1}^n(1-y_i^0)I(z_i=j)E\left(1-Y_i|{\mathit{y}}^0,\mathit{z},{\widehat{\Theta }}^{(m)}\right)\right]ln\left[p(1-\theta ){\pi }_{0j}\right]\right\}.\hfill \end{array}$$

(4)

Let ${\widehat{\mathcal{A}}}_y^{(m)}={\sum }_{j=1}^B{\widehat{\mathcal{A}}}_{yj}^{(m)}$, $y=0,1$, where

$$\begin{array}{cc}\hfill {\widehat{\mathcal{A}}}_{1j}^{(m)}& =\sum _{i=1}^{n}E\left\{Y_i{Y}_i^{0}I(Z_i=j)|Z_i=j,Y_i^0=1,{\widehat{\Theta }}^{(m)}\right\}\hfill \\ \hfill & =\sum _{i=1}^n\left[y_i^{0}I(z_i=j)\left\{\frac{p{\widehat{\theta }}^{(m)}{\widehat{\pi }}_{1j}^{(m)}}{p{\widehat{\theta }}^{(m)}{\widehat{\pi }}_{1j}^{(m)}+(1-p)(1-{\widehat{\theta }}^{(m)}){\widehat{\pi }}_{0j}^{(m)}}\right\}\right],\hfill \\ \hfill {\widehat{\mathcal{A}}}_{0j}^{(m)}& =\sum _{i=1}^{n}E\left\{Y_i(1-Y_i^0)I(Z_i=j)|Z_i=j,Y_i^0=0,{\widehat{\Theta }}^{(m)}\right\}\hfill \\ \hfill & =\sum _{i=1}^n\left[(1-y_i^0)I(z_i=j)\left\{\frac{(1-p){\widehat{\theta }}^{(m)}{\widehat{\pi }}_{1j}^{(m)}}{(1-p){\widehat{\theta }}^{(m)}{\widehat{\pi }}_{1j}^{(m)}+p(1-{\widehat{\theta }}^{(m)}){\widehat{\pi }}_{0j}^{(m)}}\right\}\right].\hfill \end{array}$$

Let

$$\begin{array}{cc}\hfill {\widehat{\mathcal{K}}}_{1j}^{(m)}& =\sum _{i=1}^{n}E\left\{(1-Y_i)Y_i^{0}I(Z_i=j)|Z_i=j,Y_i^0=1,{\widehat{\Theta }}^{(m)}\right\}\hfill \\ \hfill & =\sum _{i=1}^n\left[y_i^{0}I(z_i=j)\left\{\frac{(1-p)(1-{\widehat{\theta }}^{(m)}){\widehat{\pi }}_{0j}^{(m)}}{p{\widehat{\theta }}^{(m)}{\widehat{\pi }}_{1j}^{(m)}+(1-p)(1-{\widehat{\theta }}^{(m)}){\widehat{\pi }}_{0j}^{(m)}}\right\}\right],\hfill \\ \hfill {\widehat{\mathcal{K}}}_{0j}^{(m)}& =\sum _{i=1}^{n}E\left\{(1-Y_i)(1-Y_i^0)I(Z_i=j)|Z_i=j,Y_i^0=0,{\widehat{\Theta }}^{(m)}\right\}\hfill \\ \hfill & =\sum _{i=1}^n\left[(1-y_i^0)I(z_i=j)\left\{\frac{p(1-{\widehat{\theta }}^{(m)}){\widehat{\pi }}_{0j}^{(m)}}{(1-p){\widehat{\theta }}^{(m)}{\widehat{\pi }}_{1j}^{(m)}+p(1-{\widehat{\theta }}^{(m)}){\widehat{\pi }}_{0j}^{(m)}}\right\}\right].\hfill \end{array}$$

Note that ${\widehat{\mathcal{K}}}_{1j}^{(m)}+{\widehat{\mathcal{A}}}_{1j}^{(m)}={\sum }_{i=1}^n{y}_i^{0}I(z_i=j)$ and ${\widehat{\mathcal{K}}}_{0j}^{(m)}+{\widehat{\mathcal{A}}}_{0j}^{(m)}={\sum }_{i=1}^n(1-y_i^0)I(z_i=j)$. These imply ${\sum }_{j=1}^B({\widehat{\mathcal{K}}}_{1j}^{(m)}+{\widehat{\mathcal{A}}}_{1j}^{(m)})={\sum }_{i=1}^n{y}_i^0$ and ${\sum }_{j=1}^B({\widehat{\mathcal{K}}}_{0j}^{(m)}+{\widehat{\mathcal{A}}}_{0j}^{(m)})={\sum }_{i=1}^n(1-y_i^0)$. Hence, ${\sum }_{j=1}^B{\widehat{\mathcal{K}}}_{1j}^{(m)}={\sum }_{i=1}^n{y}_i^0-{\widehat{\mathcal{A}}}_1^{(m)}$ and ${\sum }_{j=1}^B{\widehat{\mathcal{K}}}_{0j}^{(m)}=n-{\sum }_{i=1}^n{y}_i^0-{\widehat{\mathcal{A}}}_0^{(m)}$. Based on these results, the Q-function in (4) can be re-written as follows:

$$\begin{array}{cc}\hfill Q\left(\Theta |{\widehat{\Theta }}^{(m)}\right)& =\sum _{j=1}^B\left\{{\widehat{\mathcal{A}}}_{1j}^{(m)}ln(p\theta {\pi }_{1j})+{\widehat{\mathcal{A}}}_{0j}^{(m)}ln\left[(1-p)\theta {\pi }_{1j}\right]\right.\hfill \\ \hfill & +\left.{\widehat{\mathcal{K}}}_{1j}^{(m)}ln\left[(1-p)(1-\theta )\theta {\pi }_{0j}\right]+{\widehat{\mathcal{K}}}_{0j}^{(m)}ln\left[p(1-\theta ){\pi }_{0j}\right]\right\}\hfill \\ \hfill & =\sum _{j=1}^B\left\{{\widehat{\mathcal{A}}}_{1j}^{(m)}\left[lnp+ln\theta +ln{\pi }_{1j}\right]+{\widehat{\mathcal{A}}}_{0j}^{(m)}\left[ln(1-p)+ln\theta +ln{\pi }_{1j}\right]\right.\hfill \\ \hfill & +\left.{\widehat{\mathcal{K}}}_{1j}^{(m)}\left[ln(1-p)+ln(1-\theta )+ln{\pi }_{0j}\right]+{\widehat{\mathcal{K}}}_{0j}^{(m)}\left[lnp+ln(1-\theta )+ln{\pi }_{0j}\right]\right\}\hfill \\ \hfill & =\sum _{j=1}^B\left\{({\widehat{\mathcal{A}}}_{1j}^{(m)}+{\widehat{\mathcal{K}}}_{0j}^{(m)})lnp+({\widehat{\mathcal{A}}}_{0j}^{(m)}+{\widehat{\mathcal{K}}}_{1j}^{(m)})ln(1-p)\right\}\hfill \\ \hfill & +\sum _{j=1}^B\left\{({\widehat{\mathcal{A}}}_{1j}^{(m)}+{\widehat{\mathcal{A}}}_{0j}^{(m)})ln\theta +({\widehat{\mathcal{K}}}_{1j}^{(m)}+{\widehat{\mathcal{K}}}_{0j}^{(m)})ln(1-\theta )\right\}\hfill \\ \hfill & +\sum _{j=1}^B\left\{({\widehat{\mathcal{A}}}_{1j}^{(m)}+{\widehat{\mathcal{A}}}_{0j}^{(m)})ln{\pi }_{1j}+({\widehat{\mathcal{K}}}_{1j}^{(m)}+{\widehat{\mathcal{K}}}_{0j}^{(m)})ln{\pi }_{0j}\right\}\hfill \\ \hfill & =\sum _{j=1}^B\left\{({\widehat{\mathcal{A}}}_{1j}^{(m)}+{\widehat{\mathcal{K}}}_{0j}^{(m)})lnp+({\widehat{\mathcal{A}}}_{0j}^{(m)}+{\widehat{\mathcal{K}}}_{1j}^{(m)})ln(1-p)\right\}\hfill \\ \hfill & +\left\{({\widehat{\mathcal{A}}}_1^{(m)}+{\widehat{\mathcal{A}}}_0^{(m)})ln\theta +(n-{\widehat{\mathcal{A}}}_1^{(m)}+{\widehat{\mathcal{A}}}_0^{(m)})ln(1-\theta )\right\}\hfill \\ \hfill & +\sum _{j=1}^{B-1}\left\{({\widehat{\mathcal{A}}}_{1j}^{(m)}+{\widehat{\mathcal{A}}}_{0j}^{(m)})ln{\pi }_{1j}+({\widehat{\mathcal{K}}}_{1j}^{(m)}+{\widehat{\mathcal{K}}}_{0j}^{(m)})ln{\pi }_{0j}\right\}\hfill \\ \hfill & +({\widehat{\mathcal{A}}}_{1B}^{(m)}+{\widehat{\mathcal{A}}}_{0B}^{(m)})ln\left(1-\sum _{j=1}^{B-1}{\pi }_{1j}\right)+({\widehat{\mathcal{K}}}_{1B}^{(m)}+{\widehat{\mathcal{K}}}_{0B}^{(m)})ln\left(1-\sum _{j=1}^{B-1}{\pi }_{0j}\right).\hfill \end{array}$$

(5)

M-step: The M-step is to maximize $Q(\Theta |{\widehat{\Theta }}^{(m)})$ in (5) with respect to $\Theta$ given ${\widehat{\Theta }}^{(m)}$ to update the estimate of $\Theta$ denoted by ${\widehat{\Theta }}^{(m+1)}=({\widehat{\theta }}^{(m+1)},{\widehat{\mathbf{\pi }}}_0^{(m+1)},{\widehat{\mathbf{\pi }}}_1^{(m+1)})$ by solving the equations $\partial Q(\Theta |{\widehat{\Theta }}^{(m)})/\partial \Theta =\mathbf{0}$, which are expressed as follows:

$$\begin{array}{c}\hfill {\widehat{\theta }}^{(m+1)}={\displaystyle \frac{{\widehat{\mathcal{A}}}_1^{(m)}+{\widehat{\mathcal{A}}}_0^{(m)}}{n}},{\widehat{\pi }}_{1j}^{(m+1)}={\displaystyle \frac{{\widehat{\mathcal{A}}}_{1j}^{(m)}+{\widehat{\mathcal{A}}}_{0j}^{(m)}}{{\widehat{\mathcal{A}}}_1^{(m)}+{\widehat{\mathcal{A}}}_0^{(m)}}},{\widehat{\pi }}_{0j}^{(m+1)}={\displaystyle \frac{{\widehat{\mathcal{K}}}_{1j}^{(m)}+{\widehat{\mathcal{K}}}_{0j}^{(m)}}{n-({\widehat{\mathcal{A}}}_1^{(m)}+{\widehat{\mathcal{A}}}_0^{(m)})}}.\end{array}$$

Iterate the E-step and M-step until $|{\widehat{\Theta }}^{(m+1)}-{\widehat{\Theta }}^{(m)}|<\epsilon =5\times {10}^{-5}$ in this study.

Under the null hypothesis $H_0:{\pi }_{1z}={\pi }_{0z}={\pi }_{·z}$, $z=1,2,\dots ,B$, the model in (1), called an alternative/full model, is reduced to the following null model:

$$\begin{array}{c}P(Y_i=y)={\theta }^y{(1-\theta )}^{1-y},\hfill \\ P(Z_i=z|Y_i=y)={\pi }_{·z},z=1,2,\dots ,B,y=0,1.\hfill \end{array}$$

(6)

Based on the above null model, $P(Z_i=z,Y_i=y)$ can then be expressed as

$$\begin{array}{c}\hfill P(Z_i=z,Y_i=y)={\theta }^y{(1-\theta )}^{1-y}{\pi }_{·z},z=1,2,\dots ,B,y=0,1.\end{array}$$

(7)

$P(Y_i=1|Z_i=j,Y_i^0=y^0)$ can be given by

$$\begin{array}{cc}\hfill P(Y_i=1|Z_i=j,Y_i^0=1)& ={\displaystyle \frac{p\theta }{p\theta +(1-p)(1-\theta )}},\hfill \\ \hfill P(Y_i=1|Z_i=j,Y_i^0=0)& ={\displaystyle \frac{(1-p)\theta }{(1-p)\theta +p(1-\theta )}},\hfill \\ \hfill P(Y_i=0|Z_i=j,Y_i^0=1)& ={\displaystyle \frac{(1-p)(1-\theta )}{p\theta +(1-p)(1-\theta )}},\hfill \\ \hfill P(Y_i=0|Z_i=j,Y_i^0=0)& ={\displaystyle \frac{p(1-\theta )}{(1-p)\theta +p(1-\theta )}}.\hfill \end{array}$$

Based on the arguments similar to those in EM algorithm 1, one can obtain an estimate of $\theta$ at the $(m+1)$th iteration of the EM algorithm as follows:

$$\begin{array}{c}\hfill {\tilde{\theta }}^{(m+1)}=\frac{{\tilde{A}}_1^{(m)}+{\tilde{A}}_0^{(m)}}{n},\end{array}$$

where

$$\begin{array}{cc}\hfill {\tilde{A}}_1^{(m)}& =\sum _{i=1}^n\left[I(y_i^0=1)\left\{\frac{p{\tilde{\theta }}^{(m)}}{p{\tilde{\theta }}^{(m)}+(1-p)(1-{\tilde{\theta }}^{(m)})}\right\}\right],\hfill \\ \hfill {\tilde{A}}_0^{(m)}& =\sum _{i=1}^n\left[I(y_i^0=0)\left\{\frac{(1-p){\tilde{\theta }}^{(m)}}{(1-p){\tilde{\theta }}^{(m)}+p(1-{\tilde{\theta }}^{(m)})}\right\}\right].\hfill \end{array}$$

Note that the ML estimate of ${\pi }_{·z}$ in (6) and (7) is ${\sum }_{i=1}^{n}I(z_i=z)/n$, $z=1,2,\dots ,B$.

2.2.3. Likelihood Ratio Test (LRT)

The LRT of Neyman and Pearson [49] is utilized to determine which model is appropriate for the association of the sensitive feature with an observed categorical random variable by testing the following hypotheses:

$$\begin{array}{cc}\hfill & H_0:{\pi }_{1j}={\pi }_{0j}={\pi }_{·j},\mathrm{for}\mathrm{all}j=1,2,\dots ,B,\hfill \\ \hfill & H_1:{\pi }_{1j}\ne {\pi }_{0j},\mathrm{for}\mathrm{some}j=1,2,\dots ,B.\hfill \end{array}$$

(8)

Let ${\mathbf{\pi }}_{·}=({\pi }_{·1},{\pi }_{·2},\dots ,{\pi }_{·B})$. Under $H_0$, the observed-data likelihood function of $(\theta ,{\mathbf{\pi }}_{·})$ is given by

$$\begin{array}{cc}\hfill & L_{obs}(\theta ,{\mathbf{\pi }}_{·}|{\mathit{y}}_0,\mathit{z})\hfill \\ \hfill =& \prod _{i=1}^n\left\{{\left[p\theta {\pi }_{·z_i}+(1-p)(1-\theta ){\pi }_{·z_i}\right]}^{I(y_i^0=1)}{\left[p(1-\theta ){\pi }_{·z_i}+(1-p)\theta {\pi }_{·z_i}\right]}^{I(y_i^0=0)}\right\}\hfill \\ \hfill =& \prod _{i=1}^n\prod _{j=1}^B{\left\{{\left[p\theta {\pi }_{·j}+(1-p)(1-\theta ){\pi }_{·j}\right]}^{I(y_i^0=1)}{\left[p(1-\theta ){\pi }_{·j}+(1-p)\theta {\pi }_{·j}\right]}^{I(y_i^0=0)}\right\}}^{I(z_i=j)}\hfill \\ \hfill =& \left\{\prod _{i=1}^n{\left[p\theta +(1-p)(1-\theta )\right]}^{I(y_i^0=1)}{\left[p(1-\theta )+(1-p)\theta \right]}^{I(y_i^0=0)}\right\}\left\{\prod _{i=1}^n\prod _{j=1}^B{\left[{\pi }_{·j}\right]}^{I(z_i=j)}\right\}\hfill \\ \hfill =& {\left[p\theta +(1-p)(1-\theta )\right]}^{{\sum }_{i=1}^{n}I(y_i^0=1)}{\left[p(1-\theta )+(1-p)\theta \right]}^{{\sum }_{i=1}^{n}I(y_i^0=0)}\prod _{j=1}^B{\left[{\pi }_{·j}\right]}^{{\sum }_{i=1}^{n}I(z_i=j)}.\hfill \end{array}$$

Therefore, one can use the following LRT statistic

$$\begin{array}{c}\hfill \Lambda =-2log\left\{\frac{{sup}_{\theta ,{\mathbf{\pi }}_{·}}{L}_{obs}(\theta ,{\mathbf{\pi }}_{·}|{\mathit{y}}_0,\mathit{z})}{{sup}_{\theta ,{\mathbf{\pi }}_0,{\mathbf{\pi }}_1}{L}_{obs}(\theta ,{\mathbf{\pi }}_0,{\mathbf{\pi }}_1|{\mathit{y}}_0,\mathit{z})}\right\}\stackrel{d}{⟶}{\chi }_{B-1}^2\end{array}$$

(9)

to determine whether to reject $H_0$ in (8).

2.3. A Latent Class Model Incorporating Warner’s RRT

2.3.1. Model

Let ${\left\{{\tilde{Z}}_{si}\right\}}_{s=1}^k$ be k observed ordinal categorical random variables each with values $1,2,\dots ,B$ for respondent i, where ${\tilde{Z}}_{si}$ denotes the respondent’s response to the sth selected question or statement by the LRT, which is significantly associated with the sensitive attribute at a significance level $\alpha$. Define ${\tilde{K}}_{si}=1$ if ${\tilde{Z}}_{si}\in \{1,2,\dots ,q^{*}\}$ and ${\tilde{K}}_{si}=0$ if ${\tilde{Z}}_{si}\in \{q^{*}+1,q^{*}+2,\dots ,B\}$, where $q^{*}$ is smaller than B. That is, ${\tilde{Z}}_{si}$ is dichotomized based on the criterion ${\tilde{Z}}_{si}\le q^{*}$. Based on a latent variable model, it is assumed that under the group to which an individual is known to belong, the corresponding observed/manifest variables are independent. Therefore, assume that ${\tilde{K}}_{1i},{\tilde{K}}_{2i},\dots ,{\tilde{K}}_{ki}$ given $Y_i$, $i=1,2,\dots ,n$, are independent. The aims are to estimate the parameters of the following LCM

$$\begin{array}{c}\hfill P(Y_i=y|Y_i^0,{\tilde{K}}_{1i},{\tilde{K}}_{2i},\dots ,{\tilde{K}}_{ki}),y=0,1,\end{array}$$

(10)

via an EM algorithm and to predict an individual in a sensitive group or not.

Define $P({\tilde{K}}_{si}=1|Y_i=y)={\alpha }_{ys}$, $y=0,1$, $s=1,2,\dots ,k$. Under the assumption that ${\tilde{K}}_{1i},{\tilde{K}}_{2i},\dots ,{\tilde{K}}_{ki}$ given $Y_i$ are independent, one can obtain the following results:

$$\begin{array}{cc}\hfill & P(Y_i^0=1,{\tilde{K}}_{1i}={\tilde{k}}_{1i},{\tilde{K}}_{2i}={\tilde{k}}_{2i},\dots ,{\tilde{K}}_{ki}={\tilde{k}}_{ki})\hfill \\ \hfill =& P(Y_i^0=1,{\tilde{K}}_{1i}={\tilde{k}}_{1i},{\tilde{K}}_{2i}={\tilde{k}}_{2i},\dots ,{\tilde{K}}_{ki}={\tilde{k}}_{ki}|Y_i=1)P(Y_i=1)\hfill \\ \hfill & +P(Y_i^0=1,{\tilde{K}}_{1i}={\tilde{k}}_{1i},{\tilde{K}}_{2i}={\tilde{k}}_{2i},\dots ,{\tilde{K}}_{ki}={\tilde{k}}_{ki}|Y_i=0)P(Y_i=0)\hfill \\ \hfill =& p\theta \prod _{s=1}^k{\alpha }_{1s}^{{\tilde{k}}_{si}}{(1-{\alpha }_{1s})}^{1-{\tilde{k}}_{si}}+(1-p)(1-\theta )\prod _{s=1}^k{\alpha }_{0s}^{{\tilde{k}}_{si}}{(1-{\alpha }_{0s})}^{1-{\tilde{k}}_{si}}\hfill \end{array}$$

(11)

and

$$\begin{array}{cc}\hfill & P(Y_i^0=0,{\tilde{K}}_{1i}={\tilde{k}}_{1i},{\tilde{K}}_{2i}={\tilde{k}}_{2i},\dots ,{\tilde{K}}_{ki}={\tilde{k}}_{ki})\hfill \\ \hfill =& P(Y_i^0=0,{\tilde{K}}_{1i}=k_{1i},{\tilde{K}}_{2i}={\tilde{k}}_{2i},\dots ,{\tilde{K}}_{ki}={\tilde{k}}_{ki}|Y_i=1)P(Y_i=1)\hfill \\ \hfill & +P(Y_i^0=0,{\tilde{K}}_{1i}={\tilde{k}}_{1i},{\tilde{K}}_{2i}={\tilde{k}}_{2i},\dots ,{\tilde{K}}_{ki}={\tilde{k}}_{ki}|Y_i=0)P(Y_i=0)\hfill \\ \hfill =& (1-p)\theta \prod _{s=1}^k{\alpha }_{1s}^{{\tilde{k}}_{si}}{(1-{\alpha }_{1s})}^{1-{\tilde{k}}_{si}}+p(1-\theta )\prod _{s=1}^k{\alpha }_{0s}^{{\tilde{k}}_{si}}{(1-{\alpha }_{0s})}^{1-{\tilde{k}}_{si}}.\hfill \end{array}$$

(12)

Let $\tilde{\Theta }=(\theta ,{\mathbf{\alpha }}_0,{\mathbf{\alpha }}_1)$, where ${\mathbf{\alpha }}_y=({\alpha }_{y1},{\alpha }_{y2},\dots ,{\alpha }_{yk})$, $y=0,1$. The observed data likelihood function of $\tilde{\Theta }$ given $({\mathit{Y}}^0,{\tilde{\mathit{K}}}_1,{\tilde{\mathit{K}}}_2,\dots ,{\tilde{\mathit{K}}}_k)$ is

$$\begin{array}{cc}\hfill & L_{obs}(\tilde{\Theta }|{\mathit{Y}}^0,{\tilde{\mathit{K}}}_1,{\tilde{\mathit{K}}}_2,\dots ,{\tilde{\mathit{K}}}_k)\hfill \\ \hfill =& \prod _{i=1}^n{\left\{p\theta \prod _{s=1}^k{\alpha }_{1s}^{{\tilde{K}}_{si}}{(1-{\alpha }_{1s})}^{1-{\tilde{K}}_{si}}+(1-p)(1-\theta )\prod _{s=1}^k{\alpha }_{0s}^{{\tilde{K}}_{si}}{(1-{\alpha }_{0s})}^{1-{\tilde{K}}_{si}}\right\}}^{Y_i^0}\hfill \\ \hfill & \times {\left\{(1-p)\theta \prod _{s=1}^k{\alpha }_{1s}^{{\tilde{K}}_{si}}{(1-{\alpha }_{1s})}^{1-{\tilde{K}}_{si}}+p(1-\theta )\prod _{s=1}^k{\alpha }_{0s}^{{\tilde{K}}_{si}}{(1-{\alpha }_{0s})}^{1-{\tilde{K}}_{si}}\right\}}^{1-Y_i^0},\hfill \end{array}$$

(13)

where ${\mathit{Y}}^0=(Y_1^0,Y_2^0,\dots ,Y_n^0)$ and ${\tilde{\mathit{K}}}_s=({\tilde{K}}_{s1},{\tilde{K}}_{s2},\dots ,{\tilde{K}}_{sn})$, $s=1,2,\dots ,k$. The $\mathsf{R}$ function optim or nlminb can be used to obtain the ML estimate of $\tilde{\Theta }$ in (13), but to avoid encountering the problem of divergence, we provide an EM algorithm in the following section to solve this problem.

2.3.2. EM Algorithm 2

Suppose that $(Y_i,Y_i^0,{\tilde{K}}_{1i},{\tilde{K}}_{2i},\dots ,{\tilde{K}}_{ki})$, $i=1,2,\dots ,n$, are observable. One can then express $P(Y_i=y,Y_i^0=y^0,{\tilde{K}}_{1i}={\tilde{k}}_{1i},{\tilde{K}}_{2i}={\tilde{k}}_{2i},\dots ,{\tilde{K}}_{ki}={\tilde{k}}_{ki})$, $y=0,1$, $y^0=0,1$, as follows:

$$\begin{array}{cc}\hfill P(Y_i=1,Y_i^0=1,{\tilde{K}}_{1i}={\tilde{k}}_{1i},{\tilde{K}}_{2i}={\tilde{k}}_{2i},\dots ,{\tilde{K}}_{ki}={\tilde{k}}_{ki})& =p\theta \prod _{s=1}^k{\alpha }_{1s}^{{\tilde{k}}_{si}}{(1-{\alpha }_{1s})}^{1-{\tilde{k}}_{si}},\hfill \\ \hfill P(Y_i=0,Y_i^0=1,{\tilde{K}}_{1i}={\tilde{k}}_{1i},{\tilde{K}}_{2i}={\tilde{k}}_{2i},\dots ,{\tilde{K}}_{ki}={\tilde{k}}_{ki})& =(1-p)(1-\theta )\prod _{s=1}^k{\alpha }_{0s}^{{\tilde{k}}_{si}}{(1-{\alpha }_{0s})}^{1-{\tilde{k}}_{si}},\hfill \\ \hfill P(Y_i=1,Y_i^0=0,{\tilde{K}}_{1i}={\tilde{k}}_{1i},{\tilde{K}}_{2i}={\tilde{k}}_{2i},\dots ,{\tilde{K}}_{ki}={\tilde{k}}_{ki})& =(1-p)\theta \prod _{s=1}^k{\alpha }_{1s}^{{\tilde{k}}_{si}}{(1-{\alpha }_{1s})}^{1-{\tilde{k}}_{si}},\hfill \\ \hfill P(Y_i=0,Y_i^0=0,{\tilde{K}}_{1i}={\tilde{k}}_{1i},{\tilde{K}}_{2i}={\tilde{k}}_{2i},\dots ,{\tilde{K}}_{ki}={\tilde{k}}_{ki})& =p(1-\theta )\prod _{s=1}^k{\alpha }_{0s}^{{\tilde{k}}_{si}}{(1-{\alpha }_{0s})}^{1-{\tilde{k}}_{si}}.\hfill \end{array}$$

The complete-data likelihood function of $\tilde{\Theta }$ given $(\mathit{Y},{\mathit{Y}}^0,{\tilde{\mathit{K}}}_1,{\tilde{\mathit{K}}}_2,\dots ,{\tilde{\mathit{K}}}_k)$ can be written as

$$\begin{array}{cc}\hfill & L_c(\tilde{\Theta }|\mathit{Y},{\mathit{Y}}^0,{\tilde{\mathit{K}}}_1,{\tilde{\mathit{K}}}_2,\dots ,{\tilde{\mathit{K}}}_k)\hfill \\ \hfill =& \prod _{i=1}^n{\left\{p\theta \prod _{s=1}^k{\alpha }_{1s}^{{\tilde{K}}_{si}}{(1-{\alpha }_{1s})}^{1-{\tilde{K}}_{si}}\right\}}^{Y_i^0{Y}_i}{\left\{(1-p)(1-\theta )\prod _{s=1}^k{\alpha }_{0s}^{{\tilde{K}}_{si}}{(1-{\alpha }_{0s})}^{1-{\tilde{K}}_{si}}\right\}}^{Y_i^0(1-Y_i)}\hfill \\ \hfill & \times {\left\{(1-p)\theta \prod _{s=1}^k{\alpha }_{1s}^{{\tilde{K}}_{si}}{(1-{\alpha }_{1s})}^{1-{\tilde{K}}_{si}}\right\}}^{(1-Y_i^0)Y_i}{\left\{p(1-\theta )\prod _{s=1}^k{\alpha }_{0s}^{{\tilde{K}}_{si}}{(1-{\alpha }_{0s})}^{1-{\tilde{K}}_{si}}\right\}}^{(1-Y_i^0)(1-Y_i)}.\hfill \end{array}$$

With a bit of algebra, $L_c(\tilde{\Theta }|\mathit{Y},{\mathit{Y}}^0,{\tilde{\mathit{K}}}_1,{\tilde{\mathit{K}}}_2,\dots ,{\tilde{\mathit{K}}}_k)$ can be re-expressed as

$$\begin{array}{cc}\hfill & L_c(\tilde{\Theta }|\mathit{Y},{\mathit{Y}}^0,{\tilde{\mathit{K}}}_1,{\tilde{\mathit{K}}}_2,\dots ,{\tilde{\mathit{K}}}_k)\hfill \\ \hfill =& p^{{\sum }_{i=1}^n\left[Y_i^0{Y}_i+(1-Y_i^0)(1-Y_i)\right]}{(1-p)}^{{\sum }_{i=1}^n\left[Y_i^0(1-Y_i)+(1-Y_i^0)Y_i\right]}{\theta }^{{\sum }_{i=1}^n{Y}_i}{(1-\theta )}^{{\sum }_{i=1}^n(1-Y_i)}\hfill \\ \hfill & \times \prod _{s=1}^k{\alpha }_{1s}^{{\sum }_{i=1}^n{Y}_i{\tilde{K}}_{si}}{(1-{\alpha }_{1s})}^{{\sum }_{i=1}^n{Y}_i(1-{\tilde{K}}_{si})}\prod _{s=1}^k{\alpha }_{0s}^{{\sum }_{i=1}^n(1-Y_i){\tilde{K}}_{si}}{(1-{\alpha }_{0s})}^{{\sum }_{i=1}^n(1-Y_i)(1-{\tilde{K}}_{si})}.\hfill \end{array}$$

The complete-data log-likelihood of $\tilde{\Theta }$ given $(\mathit{Y},{\mathit{Y}}^0,{\tilde{\mathit{K}}}_1,{\tilde{\mathit{K}}}_2,\dots ,{\tilde{\mathit{K}}}_k)$ can then be expressed as

$$\begin{array}{cc}\hfill & {\ell }_c(\tilde{\Theta }|\mathit{Y},{\mathit{Y}}_0,{\tilde{\mathit{K}}}_1,{\tilde{\mathit{K}}}_2,\dots ,{\tilde{\mathit{K}}}_k)\hfill \\ \hfill =& \sum _{i=1}^n\left\{[Y_i^0{Y}_i+(1-Y_i^0)(1-Y_i)]lnp+[Y_i^0(1-Y_i)+(1-Y_i^0)Y_i]ln(1-p)\right\}\hfill \\ \hfill & +\sum _{i=1}^n\left\{Y_{i}ln\theta +(1-Y_i)ln(1-\theta )\right\}\hfill \\ \hfill & +\sum _{s=1}^k\sum _{i=1}^n\left\{Y_i{\tilde{K}}_{si}ln{\alpha }_{1s}+Y_i(1-{\tilde{K}}_{si})ln(1-{\alpha }_{1s})+(1-Y_i){\tilde{K}}_{si}ln{\alpha }_{0s}\right.\hfill \\ \hfill & \left.+(1-Y_i)(1-{\tilde{K}}_{si})ln(1-{\alpha }_{0s})\right\}.\hfill \end{array}$$

(14)

$Y_i$ is a latent or unobserved variable, but one can obtain $P(Y_i=y|Y_i^0=y^0,{\tilde{K}}_{1i}={\tilde{k}}_{1i},{\tilde{K}}_{2i}={\tilde{k}}_{2i},\dots ,{\tilde{K}}_{ki}={\tilde{k}}_{ki})$, $y=0,1$, $y^0=0,1$, as follows:

$$\begin{array}{cc}\hfill & P(Y_i=1|Y_i^0=1,{\tilde{K}}_{1i}={\tilde{k}}_{1i},{\tilde{K}}_{2i}={\tilde{k}}_{2i},\dots ,{\tilde{K}}_{ki}={\tilde{k}}_{ki})\hfill \\ \hfill =& \frac{{\displaystyle p\theta \prod _{s=1}^k{\alpha }_{1s}^{{\tilde{k}}_{si}}{(1-{\alpha }_{1s})}^{1-{\tilde{k}}_{si}}}}{{\displaystyle p\theta \prod _{s=1}^k{\alpha }_{1s}^{{\tilde{k}}_{si}}{(1-{\alpha }_{1s})}^{1-{\tilde{k}}_{si}}+(1-p)(1-\theta )\prod _{s=1}^k{\alpha }_{0s}^{{\tilde{k}}_{si}}{(1-{\alpha }_{0s})}^{1-{\tilde{k}}_{si}}}},\hfill \end{array}$$

$$\begin{array}{cc}\hfill & P(Y_i=0|Y_i^0=1,{\tilde{K}}_{1i}={\tilde{k}}_{1i},{\tilde{K}}_{2i}={\tilde{k}}_{2i},\dots ,{\tilde{K}}_{ki}={\tilde{k}}_{ki})\hfill \\ \hfill =& \frac{{\displaystyle (1-p)(1-\theta )\prod _{s=1}^k{\alpha }_{0s}^{{\tilde{k}}_{si}}{(1-{\alpha }_{0s})}^{1-{\tilde{k}}_{si}}}}{{\displaystyle p\theta \prod _{s=1}^k{\alpha }_{1s}^{{\tilde{k}}_{si}}{(1-{\alpha }_{1s})}^{1-{\tilde{k}}_{si}}+(1-p)(1-\theta )\prod _{s=1}^k{\alpha }_{0s}^{{\tilde{k}}_{si}}{(1-{\alpha }_{0s})}^{1-{\tilde{k}}_{si}}}},\hfill \end{array}$$

$$\begin{array}{cc}\hfill & P(Y_i=1|Y_i^0=0,{\tilde{K}}_{1i}={\tilde{k}}_{1i},{\tilde{K}}_{2i}={\tilde{k}}_{2i},\dots ,{\tilde{K}}_{ki}={\tilde{k}}_{ki})\hfill \\ \hfill =& \frac{{\displaystyle (1-p)\theta \prod _{s=1}^k{\alpha }_{1s}^{{\tilde{k}}_{si}}{(1-{\alpha }_{1s})}^{1-{\tilde{k}}_{si}}}}{{\displaystyle (1-p)\theta \prod _{s=1}^k{\alpha }_{1s}^{{\tilde{k}}_{si}}{(1-{\alpha }_{1s})}^{1-{\tilde{k}}_{si}}+p(1-\theta )\prod _{s=1}^k{\alpha }_{0s}^{{\tilde{k}}_{si}}{(1-{\alpha }_{0s})}^{1-{\tilde{k}}_{si}}}},\hfill \end{array}$$

$$\begin{array}{cc}\hfill & P(Y_i=0|Y_i^0=0,{\tilde{K}}_{1i}=k_{1i},{\tilde{K}}_{2i}={\tilde{k}}_{2i},\dots ,{\tilde{K}}_{ki}={\tilde{k}}_{ki})\hfill \\ \hfill =& \frac{{\displaystyle p(1-\theta )\prod _{s=1}^k{\alpha }_{0s}^{{\tilde{k}}_{si}}{(1-{\alpha }_{0s})}^{1-{\tilde{k}}_{si}}}}{{\displaystyle (1-p)\theta \prod _{s=1}^k{\alpha }_{1s}^{{\tilde{k}}_{si}}{(1-{\alpha }_{1s})}^{1-{\tilde{k}}_{si}}+p(1-\theta )\prod _{s=1}^k{\alpha }_{0s}^{{\tilde{k}}_{si}}{(1-{\alpha }_{0s})}^{1-{\tilde{k}}_{si}}}}.\hfill \end{array}$$

E-step: Let ${\widehat{\tilde{\Theta }}}^{(m)}=({\widehat{\theta }}^{(m)},{\widehat{\mathbf{\alpha }}}_0^{(m)},{\widehat{\mathbf{\alpha }}}_1^{(m)})$ be an estimate of $\tilde{\Theta }$ at the mth iteration, where ${\widehat{\mathbf{\alpha }}}_y^{(m)}=({\widehat{\alpha }}_{y1}^{(m)},{\widehat{\alpha }}_{y2}^{(m)},\dots ,{\widehat{\alpha }}_{yk}^{(m)})$, $y=0,1$. Let ${\widehat{\theta }}^{(0)}=max\{0.001,[{\overline{y}}^0-(1-p)]/(2p-1)\}$, where ${\overline{y}}^0=n^{-1}{\sum }_{i=1}^n{y}_i^0$, and ${\widehat{\alpha }}_{0s}^{(0)}={\widehat{\alpha }}_{1s}^{(0)}=n^{-1}{\sum }_{i=1}^n{\tilde{k}}_{si}$ be initial values. Given observed data $({\mathit{y}}^0,{\tilde{\mathit{k}}}_1,{\tilde{\mathit{k}}}_2,\dots ,{\tilde{\mathit{k}}}_k)$, by taking the expectation of ${\ell }_c(\tilde{\Theta }|\mathit{Y},{\mathit{Y}}^0,{\tilde{\mathit{K}}}_1,{\tilde{\mathit{K}}}_2,\dots ,{\tilde{\mathit{K}}}_k)$ in (14) given $({\mathit{y}}^0,{\tilde{\mathit{k}}}_1,{\tilde{\mathit{k}}}_2,\dots ,{\tilde{\mathit{k}}}_k,{\widehat{\tilde{\Theta }}}^{(m)})$, the Q-function can be given as follows:

$$\begin{array}{c}\hfill Q\left(\Theta |{\widehat{\tilde{\Theta }}}^{(m)}\right)=E\left[{\ell }_c(\tilde{\Theta }|\mathit{Y},{\mathit{Y}}^0,{\tilde{\mathit{K}}}_1,{\tilde{\mathit{K}}}_2,\dots ,{\tilde{\mathit{K}}}_k)|{\mathit{y}}^0,{\tilde{\mathit{k}}}_1,{\tilde{\mathit{k}}}_2,\dots ,{\tilde{\mathit{k}}}_k,{\widehat{\tilde{\Theta }}}^{(m)}\right].\end{array}$$

(15)

Define ${\widetilde{\tilde{\mathcal{B}}}}_{1i}^{(m)}$ and ${\widetilde{\tilde{\mathcal{B}}}}_{0i}^{(m)}$ as follows:

$$\begin{array}{cc}\hfill {\widetilde{\tilde{\mathcal{B}}}}_{1i}^{(m)}& =E\left\{Y_i|Y_i^0=1,{\tilde{K}}_{1i}={\tilde{k}}_{1i},{\tilde{K}}_{2i}={\tilde{k}}_{2i},\dots ,{\tilde{K}}_{ki}={\tilde{k}}_{ki},{\widehat{\tilde{\Theta }}}^{(m)}\right\}\hfill \\ \hfill & =\frac{{\displaystyle I(Y_i^0=1)p{\widehat{\theta }}^{(m)}\prod _{r=1}^k{\left[{\widehat{\alpha }}_{1r}^{(m)}\right]}^{{\tilde{k}}_{ri}}{[1-{\widehat{\alpha }}_{1r}^{(m)}]}^{1-{\tilde{k}}_{ri}}}}{{\displaystyle p{\widehat{\theta }}^{(m)}\prod _{r=1}^k{\left[{\widehat{\alpha }}_{1r}^{(m)}\right]}^{{\tilde{k}}_{ri}}{[1-{\widehat{\alpha }}_{1r}^{(m)}]}^{1-{\tilde{k}}_{ri}}+(1-p)(1-{\widehat{\theta }}^{(m)})\prod _{r=1}^k{\left[{\widehat{\alpha }}_{0r}^{(m)}\right]}^{{\tilde{k}}_{ri}}{[1-{\widehat{\alpha }}_{0r}^{(m)}]}^{1-{\tilde{k}}_{ri}}}},\hfill \end{array}$$

$$\begin{array}{cc}\hfill {\widetilde{\tilde{\mathcal{B}}}}_{0i}^{(m)}& =E\left\{Y_i|Y_i^0=0,{\tilde{K}}_{1i}={\tilde{k}}_{1i},{\tilde{K}}_{2i}={\tilde{k}}_{2i},\dots ,{\tilde{K}}_{ki}={\tilde{k}}_{ki},{\widehat{\tilde{\Theta }}}^{(m)}\right\}\hfill \\ \hfill & =\frac{{\displaystyle I(Y_i^0=0)(1-p){\widehat{\theta }}^{(m)}\prod _{r=1}^k{\left[{\widehat{\alpha }}_{1r}^{(m)}\right]}^{{\tilde{k}}_{ri}}{(1-{\widehat{\alpha }}_{1r}^{(m)})}^{1-{\tilde{k}}_{ri}}}}{{\displaystyle (1-p){\widehat{\theta }}^{(m)}\prod _{r=1}^k{\left[{\widehat{\alpha }}_{1r}^{(m)}\right]}^{{\tilde{k}}_{ri}}{(1-{\widehat{\alpha }}_{1r}^{(m)})}^{1-{\tilde{k}}_{ri}}+p(1-{\widehat{\theta }}^{(m)})\prod _{r=1}^k{\left[{\widehat{\alpha }}_{0r}^{(m)}\right]}^{{\tilde{k}}_{ri}}{(1-{\widehat{\alpha }}_{0r}^{(m)})}^{1-{\tilde{k}}_{ri}}}}.\hfill \end{array}$$

Let ${\widetilde{\tilde{\mathcal{A}}}}_1^{(m)}={\sum }_{i=1}^n{\widetilde{\tilde{\mathcal{B}}}}_{1i}^{(m)}$ and ${\widetilde{\tilde{\mathcal{A}}}}_0^{(m)}={\sum }_{i=1}^n{\widetilde{\tilde{\mathcal{B}}}}_{0i}^{(m)}$. Hence, one can obtain

$$\begin{array}{cc}\hfill & Q\left(\Theta |{\widehat{\tilde{\Theta }}}^{(m)}\right)\hfill \\ \hfill =& E\left[{\ell }_c\left(\tilde{\Theta }|\mathit{Y},{\mathit{Y}}^0,{\tilde{\mathit{K}}}_1,{\tilde{\mathit{K}}}_2,\dots ,{\tilde{\mathit{K}}}_k|{\mathit{y}}^0,{\tilde{\mathit{k}}}_1,{\tilde{\mathit{k}}}_2,\dots ,{\tilde{\mathit{k}}}_k,{\widehat{\tilde{\Theta }}}^{(m)}\right)\right]\hfill \\ \hfill =& \sum _{i=1}^n\left\{\left(Y_i^0{\widetilde{\tilde{\mathcal{B}}}}_{1i}^{(m)}+(1-Y_{i0})(1-{\widetilde{\tilde{\mathcal{B}}}}_{0i}^{(m)})\right)lnp+\left(Y_i^0(1-{\widetilde{\tilde{\mathcal{B}}}}_{1i}^{(m)})+(1-Y_i^0){\widetilde{\tilde{\mathcal{B}}}}_{0i}^{(m)}\right)ln(1-p)\right\}\hfill \\ \hfill & +\sum _{i=1}^n\left\{\left({\widetilde{\tilde{\mathcal{B}}}}_{1i}^{(m)}+{\widetilde{\tilde{\mathcal{B}}}}_{1i}^{(m)}\right)ln\theta +\left[1-\left({\widetilde{\tilde{\mathcal{B}}}}_{1i}^{(m)}+{\widetilde{\tilde{\mathcal{B}}}}_{0i}^{(m)}\right)\right]ln(1-\theta )\right\}\hfill \\ \hfill & +\sum _{s=1}^k\sum _{i=1}^n\left\{\left({\widetilde{\tilde{\mathcal{B}}}}_{1i}^{(m)}+{\widetilde{\tilde{\mathcal{B}}}}_{0i}^{(m)}\right){\tilde{k}}_{si}ln{\alpha }_{1s}+\left({\widetilde{\tilde{\mathcal{B}}}}_{1i}^{(m)}+{\widetilde{\tilde{\mathcal{B}}}}_{0i}^{(m)}\right)(1-{\tilde{k}}_{si})ln(1-{\alpha }_{1s})\right\}\hfill \\ \hfill & +\sum _{s=1}^k\sum _{i=1}^n\left\{\left[1-\left({\widetilde{\tilde{\mathcal{B}}}}_{1i}^{(m)}+{\widetilde{\tilde{\mathcal{B}}}}_{0i}^{(m)}\right)\right]{\tilde{k}}_{si}ln{\alpha }_{0s}+\left[1-\left({\widetilde{\tilde{\mathcal{B}}}}_{1i}^{(m)}+{\widetilde{\tilde{\mathcal{B}}}}_{0i}^{(m)}\right)\right](1-{\tilde{k}}_{si})ln(1-{\alpha }_{0s})\right\}\hfill \\ \hfill =& \left({\widetilde{\tilde{\mathcal{A}}}}_1^{(m)}+\sum _{i=1}^n(1-Y_{i0})-{\widetilde{\tilde{\mathcal{A}}}}_0^{(m)}\right)lnp+\left(\sum _{i=1}^n{Y}_i^0-{\widetilde{\tilde{\mathcal{A}}}}_1^{(m)}+{\widetilde{\tilde{\mathcal{A}}}}_0^{(m)}\right)ln(1-p)\hfill \\ \hfill & +\left({\widetilde{\tilde{\mathcal{A}}}}_1^{(m)}+{\widetilde{\tilde{\mathcal{A}}}}_0^{(m)}\right)ln\theta +\left(n-{\widetilde{\tilde{\mathcal{A}}}}_1^{(m)}-{\widetilde{\tilde{\mathcal{A}}}}_0^{(m)}\right)ln(1-\theta )\hfill \\ \hfill & +\sum _{s=1}^k\sum _{i=1}^n\left\{\left({\widetilde{\tilde{\mathcal{B}}}}_{1i}^{(m)}+{\widetilde{\tilde{\mathcal{B}}}}_{0i}^{(m)}\right){\tilde{k}}_{si}ln{\alpha }_{1s}+\left({\widetilde{\tilde{\mathcal{B}}}}_{1i}^{(m)}+{\widetilde{\tilde{\mathcal{B}}}}_{0i}^{(m)}\right)(1-{\tilde{k}}_{si})ln(1-{\alpha }_{1s})\right\}\hfill \\ \hfill & +\sum _{s=1}^k\sum _{i=1}^n\left\{\left[1-\left({\widetilde{\tilde{\mathcal{B}}}}_{1i}^{(m)}+{\widetilde{\tilde{\mathcal{B}}}}_{0i}^{(m)}\right)\right]{\tilde{k}}_{si}ln{\alpha }_{0s}\right.\hfill \\ \hfill & \left.+\left[1-\left({\widetilde{\tilde{\mathcal{B}}}}_{1i}^{(m)}+{\widetilde{\tilde{\mathcal{B}}}}_{0i}^{(m)}\right)\right](1-{\tilde{k}}_{si})ln(1-{\alpha }_{0s})\right\}.\hfill \end{array}$$

(16)

M-step: Update the estimate of $\Theta$, denoted by ${\widehat{\tilde{\Theta }}}^{(m+1)}$ by maximizing (15), given as follows:

$$\begin{array}{cc}\hfill {\widehat{\theta }}^{(m+1)}& =\frac{{\widetilde{\tilde{\mathcal{A}}}}_1^{(m)}+{\widetilde{\tilde{\mathcal{A}}}}_0^{(m)}}{n},\hfill \\ \hfill {\widehat{\alpha }}_{1s}^{(m+1)}& =\frac{{\displaystyle \sum _{i=1}^n\left({\widetilde{\tilde{\mathcal{B}}}}_{1i}^{(m)}+{\widetilde{\tilde{\mathcal{B}}}}_{0i}^{(m)}\right){\tilde{k}}_{si}}}{{\displaystyle \sum _{i=1}^n\left({\widetilde{\tilde{\mathcal{B}}}}_{1i}^{(m)}+{\widetilde{\tilde{\mathcal{B}}}}_{0i}^{(m)}\right)}},{\widehat{\alpha }}_{0s}^{(m+1)}=\frac{{\displaystyle \sum _{i=1}^n\left[1-\left({\widetilde{\tilde{\mathcal{B}}}}_{1i}^{(m)}+{\widetilde{\tilde{\mathcal{B}}}}_{0i}^{(m)}\right)\right]{\tilde{k}}_{si}}}{{\displaystyle \sum _{i=1}^n\left[1-\left({\widetilde{\tilde{\mathcal{B}}}}_{1i}^{(m)}+{\widetilde{\tilde{\mathcal{B}}}}_{0i}^{(m)}\right)\right]}},s=1,2,\dots ,k.\hfill \end{array}$$

The E-step and M-step are iterated until $|{\widehat{\tilde{\Theta }}}^{(m+1)}-{\widehat{\tilde{\Theta }}}^{(m)}|<\epsilon =5\times {10}^{-5}$ in this study.

Based on the estimate $\widehat{\tilde{\Theta }}=(\widehat{\theta },{\widehat{\mathbf{\alpha }}}_0,{\widehat{\mathbf{\alpha }}}_1)$ of $\tilde{\Theta }=(\theta ,{\mathbf{\alpha }}_0,{\mathbf{\alpha }}_1)$, we propose the first criterion of classifying whether or not the ith individual belongs to the sensitive group, which is described below. If the ith individual were with $(Y_i^0=1,{\tilde{K}}_{1i}={\tilde{k}}_{1i},{\tilde{K}}_{2i}={\tilde{k}}_{2i},\dots ,{\tilde{K}}_{ki}={\tilde{k}}_{ki})$ and

$$\begin{array}{cc}\hfill & \widehat{P}(Y_i=1|Y_i^0=1,{\tilde{K}}_{1i}={\tilde{k}}_{1i},{\tilde{K}}_{2i}={\tilde{k}}_{2i},\dots ,{\tilde{K}}_{ki}={\tilde{k}}_{ki})\hfill \\ \hfill =& \frac{{\displaystyle p\widehat{\theta }\prod _{s=1}^k{\widehat{\alpha }}_{1s}^{{\tilde{k}}_{si}}{(1-{\widehat{\alpha }}_{1s})}^{1-{\tilde{k}}_{si}}}}{{\displaystyle p\widehat{\theta }\prod _{s=1}^k{\widehat{\alpha }}_{1s}^{{\tilde{k}}_{si}}{(1-{\widehat{\alpha }}_{1s})}^{1-{\tilde{k}}_{si}}+(1-p)(1-\widehat{\theta })\prod _{s=1}^k{\widehat{\alpha }}_{0s}^{{\tilde{k}}_{si}}{(1-{\widehat{\alpha }}_{0s})}^{1-{\tilde{k}}_{si}}}}\ge 0.5,\hfill \end{array}$$

(17)

then she/he is predicted to belong to the sensitive group. If she/he were with $(Y_i^0=0,{\tilde{K}}_{1i}={\tilde{k}}_{1i},{\tilde{K}}_{2i}={\tilde{k}}_{2i},\dots ,{\tilde{K}}_{ki}={\tilde{k}}_{ki})$ and

$$\begin{array}{cc}\hfill & \widehat{P}(Y_i=1|Y_i^0=0,{\tilde{K}}_{1i}={\tilde{k}}_{1i},{\tilde{K}}_{2i}={\tilde{k}}_{2i},\dots ,{\tilde{K}}_{ki}={\tilde{k}}_{ki})\hfill \\ \hfill =& \frac{{\displaystyle (1-p)\widehat{\theta }\prod _{s=1}^k{\widehat{\alpha }}_{1s}^{{\tilde{k}}_{si}}{(1-{\widehat{\alpha }}_{1s})}^{1-{\tilde{k}}_{si}}}}{{\displaystyle (1-p)\widehat{\theta }\prod _{s=1}^k{\widehat{\alpha }}_{1s}^{{\tilde{k}}_{si}}{(1-{\widehat{\alpha }}_{1s})}^{1-{\tilde{k}}_{si}}+p(1-\widehat{\theta })\prod _{s=1}^k{\widehat{\alpha }}_{0s}^{{\tilde{k}}_{si}}{(1-{\widehat{\alpha }}_{0s})}^{1-{\tilde{k}}_{si}}}}\ge 0.5.\hfill \end{array}$$

(18)

then she/he is classified to be in the sensitive group.

3. Real Data Application

The proposed methodology is employed to demonstrate its practical applications by using the survey study data of the sexuality of 3027 freshmen (1193 females, 1792 males, 29 non-binary, 13 no response), not including international students, enrolled based on convenience sampling at Feng Chia University in Taiwan in 2016. Because 262 (8.7%) respondents had missing data, 2765 respondents are used for the purpose. Through this face-to-face survey, respondents answered a questionnaire of 67 statements or questions, including three parts: demographic questions, attitude to love, and online dating experience. This questionnaire did not collect personal information such as respondents’ name and age. In addition, data of responses to the sensitive question, question 37: “What is your sexual orientation?”, were collected by using the multichotomous RR design of Groenitz [50,51]. We collected data of responses to the two sensitive questions, question 38: “Have you ever had sex?” and question 39: “Have you ever had a one-night stand through a dating site or mobile app?”, via the design of [3] and birthday interval as a non-random device. The generalized RR design of [10] was employed to collect data of responses to the sensitive question, question 67: “Have you ever had a one-night stand through a dating site or mobile app?”, which is the same question as question 39. Therefore, the privacy of respondents was protected.

There were 27 statements in the third part about online dating experience for respondents to answer. In addition, question 39: “Have you ever had a one-night stand through a dating site or mobile app?” was designed to collect data on the sexual behavior. Because people in eastern cultures are often shy about referring to sexual behavior or reluctant to provide correct answers, the indirect question technique of [3] was used in this study. Expressly, two birthday intervals, January 1 to August 30 and August 11 to December 31, were set up as a non-random device so that interviewees based on their birthday could answer “yes” or “no” to question 39. The design detail of question 39 is as given in Table 10 of [17]. More specifically, if a respondent whose birthday was between January 1 and August 10 had ever had a one-night stand through a dating site or mobile app, answer “A” denoted “yes” to the question, and answer “B” denoted “no”. If a respondent whose birthday was between August 11 and December 31 had ever had a one-night stand through a dating site or mobile app, answer “B” denoted “yes” to the question, and answer “A” denoted “no”. In addition, statements 40–66 are direct ones about interviewees’ thoughts and evaluation on making friends on the internet, finding a lover or sex partner, and attitude toward online dating. Each statement had five response options, i.e., $1=“\mathrm{very}\mathrm{consistent}”$, $2=“\mathrm{almost}\mathrm{consistent}”$, $3=“\mathrm{fairly}\mathrm{consistent}”$, $4=“\mathrm{a}\mathrm{bit}\mathrm{consistent}”$, and $5=“\mathrm{very}\mathrm{inconsistent}”$.

As in Section 2.2, the binary variable $Y_i$ denotes whether respondent i has ever had a one-night stand through a dating site or mobile app in which $Y_i=1$ indicates yes, and $Y_i=0$ otherwise. $P(Y_i=1)=\theta$ is the population proportion of the sensitive attribute, having ever had a one-night stand through a dating site or mobile app. $p=0.6082$ is the probability of people whose birthday is between January 1 and August 10 by assuming that birthdays are uniformly distributed. The LRT statistic in (9) is utilized to identify which statements from statements 40–66 are significantly related to question 39 based on data of responses to these 27 statements and question 39 at $\alpha =0.1$.

Table 1. Results of the likelihood ratio test and the null and alternative/full models for the association between the sensitive characteristic variable from question 39: “Have you ever had a one-night stand through a dating site or mobile app?”, and each of the 27 observed ordinal categorical random variables for online dating experience corresponding to statements 40 to 66.

Parameter	$\mathit{\theta }$	${\mathit{\pi }}_{01}$	${\mathit{\pi }}_{02}$	${\mathit{\pi }}_{03}$	${\mathit{\pi }}_{04}$	${\mathit{\pi }}_{05}$	${\mathit{\pi }}_{11}$	${\mathit{\pi }}_{12}$	${\mathit{\pi }}_{13}$	${\mathit{\pi }}_{14}$	${\mathit{\pi }}_{15}$	p-Value
Q${}_{40}$ $H_1$	0.2106	0.1083	0.1775	0.2702	0.2445	0.1995	0.1538	0.1814	0.2236	0.1275	0.3136	0.6644
$H_0$	0.2084	0.1180	0.1780	0.2600	0.2200	0.2240
Q${}_{41}$ $H_1$	0.2144	0.0401	0.0989	0.2833	0.3089	0.2689	0.0691	0.1438	0.0433	0.3829	0.3609	0.1905
$H_0$	0.2084	0.0460	0.1080	0.2320	0.3250	0.2890
Q${}_{42}$ $H_1$	0.2105	0.0226	0.1326	0.2437	0.3439	0.2571	0.1007	0.0626	0.3316	0.2168	0.2883	0.3650
$H_0$	0.2084	0.0390	0.1180	0.2620	0.3170	0.2640
Q${}_{43}$ $H_1$	0.2137	0.0347	0.1494	0.3448	0.3156	0.1555	0.1822	0.2913	0.3848	0.0590	0.0827	0.0187
$H_0$	0.2084	0.0660	0.1800	0.3530	0.2610	0.1400
Q${}_{44}$ $H_1$	0.2098	0.0335	0.1344	0.2532	0.3433	0.2356	0.0636	0.0730	0.3340	0.3290	0.2004	0.8481
$H_0$	0.2084	0.0400	0.1220	0.2700	0.3400	0.2280
Q${}_{45}$ $H_1$	0.2103	0.0704	0.2526	0.3725	0.2018	0.1028	0.1606	0.1282	0.3657	0.2467	0.0989	0.0836
$H_0$	0.2084	0.0890	0.2260	0.3710	0.2110	0.1020
Q${}_{46}$ $H_1$	0.2108	0.0753	0.2644	0.3688	0.1905	0.1009	0.1625	0.2111	0.3143	0.2114	0.1008	0.8095
$H_0$	0.2084	0.0940	0.2530	0.3570	0.1950	0.1010
Q${}_{47}$ $H_1$	0.2095	0.0805	0.2249	0.3769	0.2054	0.1123	0.1002	0.2614	0.3092	0.2297	0.0995	1.0000
$H_0$	0.2084	0.0850	0.2330	0.3630	0.2100	0.1100
Q${}_{48}$ $H_1$	0.2107	0.1024	0.2693	0.3614	0.1716	0.0954	0.2465	0.2786	0.2871	0.1073	0.0805	0.5409
$H_0$	0.2084	0.1330	0.2710	0.3460	0.1580	0.0920
Q${}_{49}$ $H_1$	0.2101	0.0297	0.1155	0.2440	0.3280	0.2828	0.0502	0.0547	0.2857	0.3041	0.3054	1.0000
$H_0$	0.2084	0.0340	0.1030	0.2530	0.3230	0.2880
Q${}_{50}$ $H_1$	0.2095	0.0368	0.1083	0.2961	0.3547	0.2040	0.0838	0.1558	0.2878	0.2394	0.2331	0.7806
$H_0$	0.2084	0.0470	0.1180	0.2940	0.3310	0.2100
Q${}_{51}$ $H_1$	0.2097	0.0353	0.1485	0.3369	0.2945	0.1847	0.1015	0.1302	0.3531	0.2042	0.2111	0.7475
$H_0$	0.2084	0.0490	0.1450	0.3400	0.2760	0.1900
Q${}_{52}$ $H_1$	0.2237	0.0230	0.1126	0.2543	0.3377	0.2724	0.0656	0.0038	0.3335	0.3315	0.2656	0.1262
$H_0$	0.2084	0.0330	0.0880	0.2720	0.3360	0.2710
Q${}_{53}$ $H_1$	0.2104	0.0356	0.1396	0.2733	0.3128	0.2386	0.0555	0.0484	0.2411	0.3369	0.3181	0.7395
$H_0$	0.2084	0.0400	0.1200	0.2670	0.3180	0.2550
Q${}_{54}$ $H_1$	0.2099	0.0269	0.0940	0.2819	0.3636	0.2336	0.0540	0.1373	0.3190	0.2182	0.2715	1.0000
$H_0$	0.2084	0.0330	0.1030	0.2900	0.3330	0.2420
Q${}_{55}$ $H_1$	0.2098	0.0062	0.0234	0.1420	0.3001	0.5282	0.0782	0.0704	0.0960	0.1953	0.5600	0.0259
$H_0$	0.2084	0.0210	0.0330	0.1320	0.2780	0.5350
Q${}_{56}$ $H_1$	0.2101	0.0104	0.0236	0.0960	0.2575	0.6124	0.0315	0.0333	0.0383	0.1681	0.7288	1.0000
$H_0$	0.2084	0.0150	0.0260	0.0840	0.2390	0.6370
Q${}_{57}$ $H_1$	0.2113	0.0758	0.2109	0.3783	0.1956	0.1393	0.2031	0.2552	0.2960	0.1993	0.0464	0.3833
$H_0$	0.2084	0.1030	0.2200	0.3610	0.1960	0.1200
Q${}_{58}$ $H_1$	0.2093	0.0328	0.0998	0.2735	0.3070	0.2870	0.0316	0.1121	0.2024	0.4333	0.2205	0.8481
$H_0$	0.2084	0.0330	0.1020	0.2590	0.3330	0.2730
Q${}_{59}$ $H_1$	0.2111	0.0355	0.0856	0.2273	0.2875	0.3641	0.0044	0.0671	0.2078	0.4214	0.2993	0.6802
$H_0$	0.2084	0.0290	0.0820	0.2230	0.3160	0.3500
Q${}_{60}$ $H_1$	0.2088	0.0278	0.0812	0.2478	0.3212	0.3220	0.0435	0.1029	0.1609	0.3263	0.3664	0.9284
$H_0$	0.2084	0.0310	0.0860	0.2300	0.3220	0.3310
Q${}_{61}$ $H_1$	0.2114	0.0002	0.0260	0.1039	0.2794	0.5904	0.0624	0.0005	0.1599	0.2253	0.5519	0.0757
$H_0$	0.2084	0.0130	0.0210	0.1160	0.2680	0.5820
Q${}_{62}$ $H_1$	0.2052	0.0000	0.0279	0.0977	0.2705	0.6039	0.0687	0.0170	0.1611	0.2390	0.5143	0.3277
$H_0$	0.2084	0.0140	0.0260	0.1110	0.2640	0.5860
Q${}_{63}$ $H_1$	0.2078	0.0001	0.0239	0.0794	0.2823	0.6143	0.0745	0.0132	0.2108	0.1630	0.5384	0.9950
$H_0$	0.2084	0.0160	0.0220	0.1070	0.2580	0.5990
Q${}_{64}$ $H_1$	0.2099	0.0030	0.0297	0.0819	0.1839	0.7015	0.0939	0.0241	0.0432	0.0953	0.7434	0.1355
$H_0$	0.2084	0.0220	0.0290	0.0740	0.1650	0.7100
Q${}_{65}$ $H_1$	0.2108	0.0035	0.0139	0.0738	0.1903	0.7184	0.0674	0.0594	0.0926	0.0733	0.7074	0.9131
$H_0$	0.2084	0.0170	0.0240	0.0780	0.1660	0.7160
Q${}_{66}$ $H_1$	0.2097	0.0048	0.0304	0.0624	0.1278	0.7747	0.0976	0.0287	0.0632	0.0859	0.7245	0.1985
$H_0$	0.2084	0.0240	0.0300	0.0630	0.1190	0.7640

displays the analysis results, including estimates of the parameters of the alternative/full model in (1) and null model in (6) for the association of this sensitive attribute with response to each of these 27 statements under $H_0$ and $H_1$ in (8) and corresponding p-values of the LRTs. As seen from these results, the p-values of the LRTs corresponding to statements 43, 45, 55, and 61 are 0.0187, 0.0836, 0.0259, and 0.0757, respectively, which imply that question 39 was significantly associated with response to each of these four statements at $\alpha =0.1$. The selected four statements are listed below.

• Q43. Finding friends on the internet can improve your social circle.
• Q45. You can find people with similar interests online.
• Q55. I am a homebody, so I want to make friends online.
• Q61. I want to find my partner through online dating.

Let ${\left\{{\tilde{Z}}_{si}\right\}}_{s=1}^4$ denote the ith respondent’s response to the aforementioned four statements selected by the LRTs at $\alpha =0.1$ with full probability model $P({\tilde{Z}}_{si}=z|Y_i=y)={\pi }_{yz}^{(s)}$, $s=1,2,3,4$, $z=1,2,3,4,5$, $y=0,1$.

Table 2. (Results extracted from Table 1) Summary of significant results of the alternative/full model for the association between the sensitive characteristic variable from question 39: “Have you ever had a one-night stand through a dating site or mobile app?”, and each of the four observed ordinal categorical random variables for online dating experience corresponding to selected statements 43, 45, 55, and 61, respectively, by the likelihood ratio test at $\alpha =0.1$.

Parameter	$\mathit{\theta }$	${\mathit{\pi }}_{01}$	${\mathit{\pi }}_{02}$	${\mathit{\pi }}_{03}$	${\mathit{\pi }}_{04}$	${\mathit{\pi }}_{05}$	${\mathit{\pi }}_{11}$	${\mathit{\pi }}_{12}$	${\mathit{\pi }}_{13}$	${\mathit{\pi }}_{14}$	${\mathit{\pi }}_{15}$
Q${}_{43}$ Estimate	0.2137	0.0347	0.1494	0.3448	0.3156	0.1555	0.1822	0.2913	0.3848	0.0590	0.0827
SE	0.0370	0.0142	0.0217	0.0225	0.0197	0.0173	0.0604	0.0721	0.0737	0.0556	0.0464
Q${}_{45}$ Estimate	0.2103	0.0704	0.2526	0.3725	0.2018	0.1028	0.1606	0.1282	0.3657	0.2467	0.0989
SE	0.0388	0.0176	0.0213	0.0320	0.0266	0.0170	0.0679	0.0642	0.1134	0.0992	0.0646
Q${}_{55}$ Estimate	0.2098	0.0062	0.0234	0.1420	0.3001	0.5282	0.0782	0.0704	0.0960	0.1953	0.5600
SE	0.0427	0.0064	0.0076	0.0174	0.0260	0.0233	0.0314	0.0299	0.0570	0.0943	0.0778
Q${}_{61}$ Estimate	0.2114	0.0002	0.0260	0.1039	0.2794	0.5904	0.0624	0.0005	0.1599	0.2253	0.5519
SE	0.0399	0.0019	0.0043	0.0156	0.0268	0.0269	0.0173	0.0058	0.0723	0.0841	0.0995

presents estimates of ${\theta }^{(s)}$, ${\pi }_{0z}^{(s)}$, and ${\pi }_{1z}^{(s)}$ with their estimated standard errors (SEs) for the alternative/full model for the association between the sensitive feature and each of ${\left\{{\tilde{Z}}_{si}\right\}}_{s=1}^4$, respectively. Note that the estimated SEs are obtained by the bootstrap method with 200 replications. It can be seen from Table 2 that the estimated SEs of ${\widehat{\theta }}^{(s)}$ are around 0.04, and the estimated SEs of ${\widehat{\pi }}_{0z}^{(s)}$ are, respectively, smaller than those of ${\widehat{\pi }}_{1z}^{(s)}$ mainly because there is more sample information to estimate ${\pi }_{0z}^{(s)}$.

Based on ${\widehat{\theta }}^{(s)}$, ${\widehat{\pi }}_{0z}^{(s)}$, and ${\widehat{\pi }}_{1z}^{(s)}$, one can estimate $P(Y_i=1|{\tilde{Z}}_{si}=z,Y_i^0=y^0)$ as follows:

$$\begin{array}{c}\hfill \widehat{P}(Y_i=1|{\tilde{Z}}_{si}=z,Y_i^0=y^0)=\left\{\begin{array}{cc}{\displaystyle \frac{p{\widehat{\theta }}^{(s)}{\widehat{\pi }}_{1z}^{(s)}}{p{\widehat{\theta }}^{(s)}{\widehat{\pi }}_{1z}^{(s)}+(1-p)(1-{\widehat{\theta }}^{(s)}){\widehat{\pi }}_{0z}^{(s)}}}\hfill & \mathrm{as}{y}^0=1,\hfill \\ {\displaystyle \frac{(1-p){\widehat{\theta }}^{(s)}{\widehat{\pi }}_{1z}^{(s)}}{(1-p){\widehat{\theta }}^{(s)}{\widehat{\pi }}_{1z}^{(s)}+p(1-{\widehat{\theta }}^{(s)}){\widehat{\pi }}_{0z}^{(s)}}}\hfill & \mathrm{as}{y}^0=0.\hfill \end{array}\right.\end{array}$$

According to the estimated posterior probability of $Y_i=1$ given ${\tilde{Z}}_{si}=z$ and $Y_i^0=y^0$, define ${\mathcal{C}}_{si}$

$$\begin{array}{c}\hfill {\mathcal{C}}_{si}=\left\{\begin{array}{cc}1\hfill & \mathrm{when}\widehat{P}(Y_i=1|{\tilde{Z}}_{si}=z,Y_i^0=y^0)\ge 0.5,\hfill \\ 0\hfill & \mathrm{when}\widehat{P}(Y_i=1|{\tilde{Z}}_{si}=z,Y_i^0=y^0)<0.5,\hfill \end{array}\right.\end{array}$$

(19)

as a conditional classifier to identify whether the ith respondent has the sensitive attribute. As a second classification criterion, the proposed classification criterion is to predict the ith respondent in the sensitive group, i.e., she/he has ever had a one-night stand through a dating website or mobile app, if any of the ${\mathcal{C}}_{si}$, $s=1,2,3,4$, is 1, i.e., if

$$\begin{array}{c}\hfill \sum _{s=1}^4{\mathcal{C}}_{si}\ge 1.\end{array}$$

(20)

By applying the proposed methodology in Section 2.3 to the real data set, we define ${\tilde{K}}_{si}=1$ if ${\tilde{Z}}_{si}\in \{1,2\}$ and ${\tilde{K}}_{si}=0$ if ${\tilde{Z}}_{si}\in \{3,4,5\}$. ${\left\{{\tilde{K}}_{si}\right\}}_{s=1}^4$ and $Y_i$ are observed binary and latent variables in the LCM, respectively. Then, $P({\tilde{K}}_{si}=1|Y_i=0)={\alpha }_{0s}$ and $P({\tilde{K}}_{si}=1|Y_i=1)={\alpha }_{1s}$, $s=1,2,3,4$.

Table 3. Estimates of parameters of the latent class model for the association between the sensitive characteristic variable from question 39: “Have you ever had a one-night stand through a dating site or mobile app?”, and the four observed dichotomized random variables for online dating experience corresponding to selected statements 43, 45, 55, and 61 by the likelihood ratio test at $\alpha =0.1$.

	$\mathit{\theta }$	${\mathit{\alpha }}_{01}$	${\mathit{\alpha }}_{02}$	${\mathit{\alpha }}_{03}$	${\mathit{\alpha }}_{04}$	${\mathit{\alpha }}_{11}$	${\mathit{\alpha }}_{12}$	${\mathit{\alpha }}_{13}$	${\mathit{\alpha }}_{14}$
Estimate	0.2201	0.0594	0.1602	0.0145	0.0052	0.9068	0.8668	0.1966	0.1361
SE	0.0127	0.0097	0.0095	0.0033	0.0019	0.0262	0.0257	0.0166	0.0142

displays the estimates of $\theta$, ${\alpha }_{0s}$, and ${\alpha }_{1s}$ in the proposed LCM incorporating the RRT of [3] with their estimated SEs. The SEs are estimated by using the bootstrap method with 200 replications. The estimate of $\theta$ is 0.2201, and its estimated SE is 0.0127. This estimate of $\theta$ is reasonable because this estimate is still within each of the 95% confidence intervals (CIs) of ${\widehat{\theta }}^{(s)}$, $({\widehat{\theta }}^{(s)}-1.96{\widehat{\mathrm{SE}}}_{{\widehat{\theta }}^{(s)}},{\widehat{\theta }}^{(s)}+1.96{\widehat{\mathrm{SE}}}_{{\widehat{\theta }}^{(s)}})$, for the selected statements 43, 45, 55, and 61, which are $(0.1412,0.2862)$, $(0.1343,0.2863)$, $(0.1261,0.2935)$, and $(0.1332,0.2896)$, respectively, where ${\widehat{\theta }}^{(s)}$ and ${\widehat{\mathrm{SE}}}_{{\widehat{\theta }}^{(s)}}$, $s=1,2,3,4$, are given in Table 2. On the other hand, if using the estimate of $\theta$, 0.2201, and its estimated SE, 0.0127, to construct a 95% CI, the 95% CI $(0.1952,0.2450)$ also contains each of these four estimates of $\theta$, $0.2137,0.2103,0.2098$, and $0.2201$, in Table 2.

Table 4. Comparison of results of predicting respondents to have ever had a one-night stand through a dating site or mobile app based on classification criterion 1 in (17) and (18) and criterion 2 in (20).

	Classification Criterion 1
Classification Criterion 2	Having Ever Had a One-Night Stand through a Dating Site or Mobile App	Having Never Had a One-Night Stand through a Dating Site or Mobile App	Total
Having ever had a one-night stand	212	34	246
through a dating site or mobile app
Having never had a one-night stand	341	2178	2519
through a dating site or mobile app
Total	553	2212	2765

shows the classifications of 2765 respondents based on the two proposed classification criteria in (17), (18) and (20), respectively. Based on classification criterion 1, the estimated proportion of respondents who “have ever had a one-night stand through a dating site or mobile app” is $0.2(553/2765)$, which is quite close to the estimates of $\theta$, 0.2201, given in Table 3, obtained from fitting the LCM in (10). By applying the second classification criterion, 246 respondents (with estimated proportion $0.089=246/2765$) are predicted to have the sensitive attribute, “having ever had a one-night stand through a dating site or mobile app”, while the estimates of $\theta$ given in Table 2 are around 0.21, which are obtained from fitting the alternative/full model for the association of this sensitive attribute with each of the four observed ordinal categorical variables for online dating experience corresponding to selected statements 43, 45, 55, and 61, respectively, by the LRT at $\alpha =0.1$.

4. Conclusions

A combination of the RRT of [3] and an LCM has been proposed to investigate the association between a sensitive character or latent variable and the observed binary random variables, which were obtained from dichotomizing the observed ordinal categorical variables selected by the LRT. The concept of the relationship between a sensitive attribute variable and auxiliary random variables in [19] has been extended by applying the RR design of [3] to collect sensitive characteristic information. The EM algorithm, called EM algorithm 1, has been provided to easily estimate the parameters of the null and alternative models for the association of the sensitive attribute variable with each of the auxiliary ordinal categorical random variables. The LRT has been utilized to select ordinal categorical random variables that are significantly associated with the sensitive attribute variable.

An LCM has been proposed to relate observed binary random variables to a sensitive characteristic or latent variable under the RRT of [3]. The EM algorithm, called EM algorithm 2, has been proposed to easily estimate its parameters and the population proportion of the sensitive characteristic. Two classification criteria have been conducted to predict the presence of the sensitive attribute in an individual. Practical applications of the proposed methodology have been demonstrated with the survey data from the study of the sexuality of freshmen, except international students, at Feng Chia University in Taiwan in 2016. Finally, the proposed methodology can be generalized to other RRTs to make inferences for a sensitive characteristic variable. This issue can also be solved by using the Bayesian approach.