2. Bayesian latent class imputation model with structural zeros

Daniel Manrique-Vallier and Jerome P. Reiter

Previous | Next

Suppose that we have a sample of $n$ individuals measured on $J$ categorical variables. Each individual has an associated response vector $x_i=\left(x_{i1}\mathrm{,}{x}_{i2}\mathrm{,}\dots \mathrm{,}{x}_{iJ}\right),$ whose components take values from a set of $L_j$ levels. For convenience, we label these levels using consecutive numbers, $x_{ij}\in \left\{\mathrm{1,}\dots \mathrm{,}{L}_j\right\},$ so that $x_i\in C={\displaystyle {\prod }_{j=1}^J}\left\{\mathrm{1,}\dots \mathrm{,}{L}_j\right\}$ Note that $C$ includes all combinations of the $J$ variables, including structural zeros, and that each combination $x$ can be viewed as a cell in the contingency table formed by $C$. Let $x_i=\left(x_i^{\text{obs}}\mathrm{,}{x}_i^{\text{mis}}\right),$ where $x_i^{\text{obs}}$ includes the variables with observed values and $x_i^{\text{mis}}$ includes the variables with missing values. Finally, let $S=\left\{s_1\mathrm{,}\dots ,s_C\right\},$ where $s_c\in C$ and $c=\mathrm{1,}\dots \mathrm{,}C<\left|S\right|,$ be the set of structural zero cells, i.e., $\mathrm{Pr}\left(x_i\in S\right)=0.$

2.1 Latent class models

As an initial step, we describe the Bayesian latent class model without any concerns for structural zeros and without any missing data, i.e., $x_i=x_i^{\text{obs}}.$ This model is a finite mixture of product-multinomial distributions,

$$p\left(x|\lambda \mathrm{,}\pi \right)=f^{\text{LCM}}\left(x|\lambda \mathrm{,}\pi \right)={\displaystyle \sum _{k=1}^K}{\pi }_k{\displaystyle \prod _{j=1}^J}{\lambda }_{jk}\left[x_j\right],\text{ (2}\text{.1)}$$

where $\lambda =\left({\lambda }_{jk}\left[l\right]\right),$ with all ${\lambda }_{jk}\left[l\right]>0$ and ${\displaystyle {\sum }_{l=1}^{L_j}}{\lambda }_{jk}\left[l\right]=1.$ Here, $\pi =\left({\pi }_1\mathrm{,}\dots \mathrm{,}{\pi }_K\right)$ with ${\displaystyle {\sum }_{k=1}^K}{\pi }_k=1.$ This model corresponds to the generative process,

$$x_{ij}|z_i\stackrel{\text{indep}}{{\displaystyle \sim }}{\text{Discrete}}_{1:L_j}\left({\lambda }_{j{z}_i}\left[1\right]\mathrm{,}\dots \mathrm{,}{\lambda }_{j{z}_i}\left[L_j\right]\right)\text{for}\text{all}i\text{and}j\text{ (2}\text{.2)}$$

$$z_i|\pi \stackrel{\text{iid}}{{\displaystyle \sim }}{\text{Discrete}}_{1:K}\left({\pi }_1\mathrm{,}\dots \mathrm{,}{\pi }_K\right)\text{for}\text{all}i\mathrm{.}\text{ (2}\text{.3)}$$

As notation, let $\left(X\mathrm{,}Z\right)$ be a sample of $n$ variates obtained from this process, with $X=\left(x_1\mathrm{,}\dots \mathrm{,}{x}_n\right)$ and $Z=\left(z_1\mathrm{,}\dots \mathrm{,}{z}_n\right).$ For $K$ large enough, (2.1) can represent arbitrary joint distributions for $x$ (Suppes and Zanotti 1981; Dunson and Xing 2009). And, using the conditional independence representation in (2.2) and (2.3), the model can be estimated and simulated from efficiently even for large $J.$

For prior distributions on $\pi ,$ we follow Si and Reiter (2013) and Manrique-Vallier and Reiter (forthcoming 2014). We have

$${\lambda }_{jk}[\cdot ]\stackrel{\text{indep}}{{\displaystyle \sim }}\text{Dirichlet}\left(1_{L_j}\right)\text{ (2}\text{.4)}$$

$${\pi }_k=V_k{\displaystyle \prod _{h<k}}\left(1-V_h\right)\text{ (2}\text{.5)}$$

$$V_k\stackrel{\text{iid}}{{\displaystyle \sim }}\text{Beta}\left(\mathrm{1,}\alpha \right)\text{for}k=\mathrm{1,}\dots \mathrm{,}K-1;V_K=1\text{ (2}\text{.6)}$$

$$\alpha \sim \text{Gamma}\left(\mathrm{0.25,0.25}\right)\text{ (2}\text{.7)}$$

The prior distributions in (2.4) are equivalent to uniform distributions over the support of the $J\times K$ multinomial conditional probabilities and hence represent vague prior knowledge. The prior distribution for $\pi$ in (2.5)-(2.7) is an example of a finite-dimensional stick-breaking prior distribution (Sethuraman 1994; Ishwaran and James 2001). As discussed in Dunson and Xing (2009) and Si and Reiter (2013), it typically allocates $Z$ to fewer than $K$ classes, thereby reducing computation and avoiding over-fitting. For further discussion and justification of this model as an imputation engine, see Si and Reiter (2013).

2.2 Truncated latent class models

The latent class model in (2.1) does not naturally specify cells with structural zeros a priori, because it assumes a positive probability for each cell. Thus, to represent tables with structural zeros, we need to truncate the model so that

$$f^{\text{TLCM}}\left(x|\lambda \mathrm{,}\pi \mathrm{,}S\right)\propto 1\left\{x\notin S\right\}{\displaystyle \sum _{k=1}^K}{\pi }_k{\displaystyle \prod _{j=1}^J}{\lambda }_{jk}\left[x_j\right]\mathrm{.}\text{ (2}\text{.8)}$$

As Manrique-Vallier and Reiter (forthcoming 2014) show, obtaining samples from the posterior distribution of parameters $\left(\lambda \mathrm{,}\pi \right),$ conditional on a sample $X^1=\left(x_1\mathrm{,}\dots \mathrm{,}{x}_n\right),$ can be greatly facilitated by adopting a sample augmentation strategy akin to those in Basu and Ebrahimi (2001) and O’Malley and Zaslavsky (2008). We consider $X^1$ to be the portion of variates that did not fall into the set $S$ from a larger sample, $X,$ generated directly from (2.1). Let $n_0,$ $X^0,$ and $Z^0$ be the the (unknown) sample size, response vectors, and latent class labels for the portion of $X$ that did fall into $S.$ Using a prior distribution from Meng and Zaslavsky (2002), Manrique-Vallier and Reiter (forthcoming 2014) show that if $p\left(N\right)\propto 1/N,$ where $N=n_0+n,$ the posterior distribution of $\left(\lambda \mathrm{,}\pi \right)$ under the truncated model (2.8) can be obtained by integrating the posterior distribution under the augmented sample model over $\left(n_0\mathrm{,}{X}^0\mathrm{,}{Z}^0\mathrm{,}{Z}^1\right).$

In doing so, Manrique-Vallier and Reiter (forthcoming 2014) develop a computationally efficient algorithm for dealing with large sets of structural zeros when they can be expressed as the union of sets defined by margin conditions. These are sets defined by fixing some levels of a subset of the categorical variables, for example, the set of all cells such that $\left\{x\in C:x_3=\mathrm{1,}{x}_6=3\right\}.$ Manrique-Vallier and Reiter (forthcoming 2014) introduce a vector notation to denote margin conditions, which we use here as well. Let $\mu =\left({\mu }_1\mathrm{,}{\mu }_2\mathrm{,}\dots \mathrm{,}{\mu }_J\right)$ where, for $j=\mathrm{1,}\dots \mathrm{,}J,$ we let ${\mu }_j=x_j$ whenever $x_j$ is fixed at some level and ${\mu }_j=\ast$ otherwise, where $\ast$ is special notation for a placeholder. Using this notation and assuming $J=8,$ the conditions that define the example set above ( $x_3=1$ and $x_6=3$ ) correspond to the vector $\left(\ast \mathrm{,}\ast \mathrm{,1,}\ast \mathrm{,}\ast \mathrm{,3,}\ast \mathrm{,}\ast \right).$ To avoid cluttering the notation, we use the vectors $\mu$ to represent both the margin conditions and the cells defined by those margin conditions, determined from context.

2.3 Estimation and multiple imputation

We now discuss how the model in Section 2.2 can be estimated, and subsequently converted into a multiple imputation engine, when some items are missing at random. The basic strategy is to use a Gibbs sampler. Given a completed dataset $\left(x^{\text{obs}}\mathrm{,}{x}^{\text{mis}}\right),$ we take a draw of the parameters using the algorithm from Manrique-Vallier and Reiter (forthcoming 2014). Given a draw of the parameters, we take a draw of $x^{\text{mis}}$ as described below.

Formally, the algorithm proceeds as follows. Suppose that the set of structural zeros can be defined as the union of $C$ disjoint margin conditions, $S={\cup }_{c=1}^C{\mu }_c,$ and that we use the priors for $\alpha ,$ $\lambda$ and $\pi$ defined in Section 2.1. Given $x_i=\left(x_i^{\text{obs}}\mathrm{,}{x}_i^{\text{mis}}\right)$ for $i=\mathrm{1,}\dots \mathrm{,}n,$ the algorithm of Manrique-Vallier and Reiter (forthcoming 2014) samples parameters as follows.

1. For $i=\mathrm{1,}\dots \mathrm{,}n,$ sample $z_i^1\sim {\text{Discrete}}_{1:K}\left(p_1\mathrm{,}\dots \mathrm{,}{p}_k\right)\mathrm{,}$with $p_k\propto {\pi }_k{\displaystyle {\prod }_{j=1}^J}{\lambda }_{jk}\left[x_{ij}^1\right].$

2. For $j=\mathrm{1,}\dots \mathrm{,}J$ and $k=\mathrm{1,}\dots \mathrm{,}K,$ sample ${\lambda }_{jk[\cdot ]}\sim \text{Dirichlet}\left({\xi }_{jk1}\mathrm{,}\dots \mathrm{,}{\xi }_{jk{L}_j}\right),$ with ${\xi }_{jkl}=1+{\displaystyle {\sum }_{i=1}^n}1\left\{x_{ij}^1=l\mathrm{,}{z}_i^1=k\right\}+{\displaystyle {\sum }_{i=1}^{n_0}}1\left\{x_{ij}^0=l\mathrm{,}{z}_i^0=k\right\}.$

3. For $k=\mathrm{1,}\dots \mathrm{,}K-1$ sample $V_k\sim \text{Beta}\left(1+{\nu }_k\mathrm{,}a+{\displaystyle {\sum }_{h=k+1}^K}{\nu }_k\right)$ where ${\nu }_k={\displaystyle {\sum }_{i=1}^n}1\left\{z_i^1=k\right\}+{\displaystyle {\sum }_{i=1}^{n_0}}1\left\{z_i^0=k\right\}.$ Let $V_K=1$ and make ${\pi }_k=V_k{\displaystyle {\prod }_{h<k}}\left(1-V_h\right)$ for all $k=\mathrm{1,}\dots \mathrm{,}K.$

4. For $c=\mathrm{1,}\dots \mathrm{,}C,$ compute ${\omega }_c=\mathrm{Pr}\left(x\in {\mu }_c|\lambda ,\pi \right)={\displaystyle {\sum }_{k=1}^K}{\pi }_k{\displaystyle {\prod }_{{\mu }_{cj}\ne \ast }}{\lambda }_{jk}\left[{\mu }_{cj}\right].$

5. Sample $\left(n_1\mathrm{,}\dots \mathrm{,}{n}_C\right)\sim NM\left(n\mathrm{,}{\omega }_1\mathrm{,}\dots \mathrm{,}{\omega }_C\right),$ where $NM$ is the negative multinomial distribution, and let $n_0={\displaystyle {\sum }_{c=1}^C}{n}_c.$

6. Let $\kappa \leftarrow 1.$ Repeat the following for each $c=\mathrm{1,}\dots \mathrm{,}C.$

(a) Compute the normalized vector $\left(p_1\mathrm{,}\dots \mathrm{,}{p}_K\right),$ where $p_k\propto {\pi }_k{\displaystyle {\prod }_{j:{\mu }_{cj}\ne \ast }{\lambda }_{jk}\left[{\mu }_{cj}\right]}.$

(b) Repeat the following three steps $n_c$ times:

i. Sample $z_{\kappa }^0\sim \text{Discrete}\left(p_1\mathrm{,}\dots \mathrm{,}{p}_k\right),$

ii. For $j=\mathrm{1,}\dots \mathrm{,}J$ sample
$$x_{\kappa j}^0\sim (\begin{array}{ll}{\text{Discrete}}_{1:L_j}\left({\lambda }_{j{z}_{\kappa }^0}\left[1\right]\mathrm{,}\dots \mathrm{,}{\lambda }_{j{z}_{\kappa }^0}\left[L_j\right]\right)\hfill & \text{if}{\mu }_{cj}=\ast \hfill \\ {\delta }_{{\mu }_{jc}}\hfill & \text{if}{\mu }_{cj}\ne \ast \hfill \end{array}$$
where ${\delta }_{{\mu }_{cj}}$ is a point mass distribution at ${\mu }_{cj},$

iii. Let $\kappa \leftarrow \kappa +1.$

7. Sample $\alpha \sim \text{Gamma}\left(a-1+K\mathrm{,}b-\log {\pi }_K\right).$

Having sampled parameters, we now need to take a draw of $x^{\text{mis}}.$ For $i=\mathrm{1,}\dots \mathrm{,}n,$ let $m_i=\left(m_{i1}\mathrm{,}\dots \mathrm{,}{m}_{iJ}\right)$ be a vector such that $m_{ij}=1$ if component $j$ in $x_i$ is missing and $m_{ij}=0$ otherwise. Assuming that data are missing at random, we need to sample only the components of each $x_i$ for which $m_{ij}=1,$ conditional on the components for which $m_{ij}=0.$ Thus, we add an eighth step to the algorithm.

8. For $i=\mathrm{1,}\dots \mathrm{,}n,$ sample $x_i^{\text{mis}}$ from its full conditional distribution,
$$p\left(x_i^{\text{mis}}|\dots \right)\propto 1\left\{x_i\notin S\right\}{\displaystyle \prod _{j:m_{ij}=1}}{\lambda }_{j{z}_i}\left[x_{ij}\right]\mathrm{.}\text{ (2}\text{.9)}$$

In the absence of structural zeros, the $x_{ij}$ to be imputed are conditionally independent given $z_i,$ making the imputation task a routine multinomial sampling exercise (Si and Reiter 2013). However, the structural zeros in $S$ induce dependency between the components. Thus, we cannot simply sample the components independently of one another. A naive approach is to use an acceptance-rejection scheme, sampling repeatedly from the proposal distribution $p\left(x^{\text{mis}\ast }\right)={\displaystyle {\prod }_{j:m_{ij}=1}}{\lambda }_{j{z}_i}\left[x_{ij}\right]$ until obtaining a variate such that $x^{\text{mis}\ast }\notin S.$ However, when the rejection region is large or has a high probability, this approach can be very inefficient.

Instead we suggest forming additional Gibbs sampling steps, computing the conditional distributions of all missing components so that they can be sampled individually. Let $\text{Rep}\left(x_i\mathrm{,}j\mathrm{,}l\right)$ be the vector that results from replacing component $j$ in $x_i$ by an arbitrary value $l\in \left\{\mathrm{1,2},\dots ,L_j\right\}.$ The full conditional distribution of missing component $j$ of $x_i$ (when $m_{ij}=1$) is $p\left(x_{ij}|\dots \right)\propto 1\{\text{Rep}\left(x_i\mathrm{,}j\mathrm{,}{x}_{ij}\right)\notin S\}{\lambda }_{j{z}_i}\left[x_{ij}\right].$ Thus, we replace step 8 in the algorithm with

8’. For each $\left(i\mathrm{,}j\right)\in \left\{\left(i\mathrm{,}j\right):m_{ij}=1\right\},$ sample $x_{ij}\sim {\text{Discrete}}_{1:L_j}\left(p_1\mathrm{,}\dots \mathrm{,}{p}_{L_j}\right),$ where $p_l\propto {\lambda }_{j{z}_i}\left[l\right]1\left\{\text{Rep}\left(x_i\mathrm{,}j\mathrm{,}l\right)\notin S\right\}.$

The definition of $p_l$ implies trimming the support of the full conditional distribution of$x_{ij}$ from $\left\{\mathrm{1,}\dots \mathrm{,}{L}_j\right\}$ to only values that avoid $x_i\in S,$ given current values of $\left\{x_{i{j}^{\prime }}:\text{all}{j}^{\prime }\ne j\right\}.$

To obtain $M$ completed datasets for use in multiple imputation, analysts select $M$ of the sampled $x^{\text{mis}}$ after convergence of the Gibbs sampler. These datasets should be spaced sufficiently so as to be approximately independent (given $x^{\text{obs}}$). This involves thinning the MCMC samples so that the autocorrelations among parameters are close to zero.

Previous | Next