An accelerated EM algorithm for mixture models with uncertainty for rating data

Abstract

The paper is framed within the literature around Louis’ identity for the observed information matrix in incomplete data problems, with a focus on the implied acceleration of maximum likelihood estimation for mixture models. The goal is twofold: to obtain direct expressions for standard errors of parameters from the EM algorithm and to reduce the computational burden of the estimation procedure for a class of mixture models with uncertainty for rating variables. This achievement fosters the feasibility of best-subset variable selection, which is an advisable strategy to identify response patterns from regression models for all Mixtures of Experts systems. The discussion is supported by simulation experiments and a real case study.

Appendices

Appendix: the EM algorithm for CUB models

Given the notation set in Sects. 2 and 4, for a sample $\varvec{R} = (R_1,\ldots ,R_n)^{\prime }$ of ordinal scores to be collected on a scale with m categories, consider the full cub specification with covariates given in (1)–(2). Then, $\varvec{R}$ denotes the so-called incomplete data; let $\varvec{X} = (\varvec{R}^{\prime }, \varvec{Z}^{\prime })^{\prime }$ be the complete data, with missing data $\varvec{Z} = (Z_{1},\ldots , Z_n)^{\prime }$ given by:

$$\begin{aligned} Z_{i}= {\left\{ \begin{array}{ll} 1 &{}\quad \text {if the i-th observation is drawn from the feeling component} \\ 0 &{}\quad \text {otherwise}. \end{array}\right. } \end{aligned}$$

To be more specific, one should set $\varvec{Z}_{1} = \varvec{Z}$ and $\varvec{Z}_{2} = 1 - \varvec{Z}_1$ for the uncertainty component. Then, with obvious notation, consider the complete log-likelihood:

$$\begin{aligned} l_c(\varvec{\theta }; \varvec{R}, \varvec{Z}) = \sum _{i=1}^n Z_{1i}\, \log \big (\pi _i\;b_{R_i}(\xi _i) \big ) + \sum _{i=1}^n (1-Z_{1i})\, \log \big ((1-\pi _i)\dfrac{1}{m} \big ). \end{aligned}$$

(23)

At the k-th iteration and for realization $\varvec{r} = (r_1,\ldots ,r_n)^{\prime }$ of $\varvec{R}$, the procedure first computes the posterior probabilities that each observation is drawn from each component: In particular:

$$\begin{aligned} \tau _{1i}^{(k)}= & {} \dfrac{\pi _i^{(k)} b_{r_i}(\xi _i^{(k)})}{Pr(R_i=r_i \vert \varvec{\theta }^{(k)}, \varvec{y}_i, \varvec{w}_i)}, \qquad \end{aligned}$$

(24)

$$\begin{aligned} \tau _{2\,i}^{(k)}= & {} 1- \tau _{1\,i}^{(k)} = \dfrac{1}{m}\dfrac{1-\pi _i^{(k)}}{Pr(R_i=r_i \vert \varvec{\theta }^{(k)}, \varvec{y}_i, \varvec{w}_i)} \end{aligned}$$

(25)

where one sets:

$$\begin{aligned} {\text {logit}}(\pi _i^{(k)}) = \bar{\varvec{y}_i}\, \varvec{\beta }^{(k)} \qquad {\text {logit}}(\xi _i^{(k)}) = \bar{\varvec{w}_i}\, \varvec{\gamma }^{(k)} \;\qquad i=1,\ldots ,n. \end{aligned}$$

Thus, at the k-th step, the conditional expectation of the complete log-likelihood (23) to be maximized over $\varvec{\theta }$ is given by:

$$\begin{aligned} Q(\varvec{\theta }; \varvec{\theta }^{(k)})= & {} {\mathbb {E}}_{\varvec{\theta }^{(k)}}[l_c(\varvec{\theta }; \varvec{R}, \varvec{Z})| \varvec{R} = \varvec{r}] \\= & {} \sum _{i=1}^{n} \tau _{1i}^{(k)} \log (\pi _i(\varvec{\beta })) + \sum _{i=1}^{n}(1- \tau _{1i}^{(k)}) \log (1 - \pi _i(\varvec{\beta })) \\&\qquad + \sum _{i=1}^{n} \tau _{1i}^{(k)} \log (b_{r_i}(\xi _i(\varvec{\gamma }))) \; + \;\sum _{i=1}^{n} (1- \tau _{1i}^{(k)}) \log \big (\dfrac{1}{m}\big ) \end{aligned}$$

yielding to the updated estimate $\varvec{\theta }^{(k+1)}$. These steps are iterated until convergence is attained within a certain numerical tolerance. The Nelder–Mead algorithm is considered for the optimization steps in $\varvec{\theta } = (\varvec{\beta },\varvec{\gamma })^{\prime }$.

Appendix: Louis’ identity for CUB models

In full generality, consider cub models specification with $p\ge 0$, $q\ge 0$ covariates for uncertainty and feeling parameters, respectively. Notice that, since ${\text {logit}}(\pi _i) = \bar{\varvec{y}}_i \cdot \varvec{\beta }$ and ${\text {logit}}(\xi _i) = \bar{\varvec{w}}_i \cdot \varvec{\gamma }$, then

$$\begin{aligned} \dfrac{\partial \pi _i}{\beta _j} = {\bar{y}}_{ij}\,\pi _i\,(1-\pi _i), \qquad \dfrac{\partial \xi _i}{\gamma _j} = {\bar{w}}_{ij}\,\xi _i\,(1-\xi _i). \end{aligned}$$

For the complete information matrix (17), the first derivatives with respect to $\beta _j$ and $\gamma _l$ of the complete log-likelihood (23) are given by:

$$\begin{aligned} \dfrac{\partial \ell _c(\varvec{\theta })}{\partial \beta _j} = \sum _{i=1}^n {\bar{y}}_{ij} \big (Z_{i1} - \pi _i\big ), \qquad \dfrac{\partial \ell _c(\varvec{\theta })}{\partial \gamma _l} = \sum _{i=1}^n {\bar{w}}_{il} \,Z_{i1}\big ( m - R_i - \xi _i\,(m-1) \big )\quad \end{aligned}$$

(26)

from which it follows:

$$\begin{aligned}&\dfrac{\partial ^2 \ell _c(\varvec{\theta })}{\partial \beta _j\, \partial \beta _k} = - \sum _{i=1}^n {\bar{y}}_{ij}\,{\bar{y}}_{ik} \,\pi _i\,(1-\pi _i), \\&\quad \dfrac{\partial ^2 \ell _c(\varvec{\theta })}{\partial \gamma _h\partial \gamma _l} = - (m-1) \sum _{i=1}^n {\bar{w}}_{il} {\bar{w}}_{ih}\,Z_{i1}\,\xi _i (1-\xi _i) \end{aligned}$$

Then, taking the conditional expectation given $\varvec{R} = \varvec{r}$ of the negative of these second-order derivatives, the block-wise definition in (17) follows. Then, the matrix specification (19) in Louis’ identity follows straightforwardly since ${\mathbb {E}}[Z_{i1}|\varvec{R}=\varvec{r}] = \tau _i$.

Starting from the complete score vector (26), matrix ${\mathcal {V}}_c$ in (18) can be obtained as follows:

$$\begin{aligned} {\mathcal {V}}_c[j,l]&= {\mathbb {E}}\big [ \big ( \sum _{i=1}^n {\bar{y}}_{ij}(Z_{i1} - \pi _i)\big )\big (\sum _{t=1}^n {\bar{y}}_{tl}(Z_{t1} - \pi _t)\big ) \big | \varvec{R} = \varvec{r} \big ] \\&= \sum _{i=1}^n \sum _{t=1}^n {\bar{y}}_{ij}{\bar{y}}_{tl}Cov(Z_{i1}-\pi _i,Z_{t1} - \pi _t\big | \varvec{R} = \varvec{r} \big ) + \\&\quad + {\mathbb {E}}\big [\sum _{i=1}^n {\bar{y}}_{ij}(Z_{i1} - \pi _i)\big ]\,{\mathbb {E}}\big [\sum _{t=1}^n {\bar{y}}_{tl}(Z_{t1} - \pi _t)\big ] \\&= \sum _{i=1}^n {\bar{y}}_{ij}{\bar{y}}_{il} \tau _{i1}(1-\tau _{i1}) + \big (\sum _{i=1}^n {\bar{y}}_{ij}(\tau _{i1}-\pi _i)\big ) \big (\sum _{i=1}^n {\bar{y}}_{il}(\tau _{i1}-\pi _i)\big ) \\&= {\varvec{Y}}_{\varvec{\tau }}[,j]\cdot \varvec{Y}[,l] + {\mathcal {V}}[j,l] \end{aligned}$$

according to the notation introduced in Sect. 4 (notice that in the second to last row of the above identity, one uses the fact that the $Z_{i1}$’s are independent Bernoulli distributed, given $\varvec{R}=\varvec{r}$, with probability parameter given by $\tau _{1i}$. Similar steps can be easily pursued for the other blocks of the matrix (18).

Finally, the score vector for the incomplete data problem is obtained by taking the first partial derivatives in (3) with respect to $\beta _j$, $j=0,\ldots ,p$ and $\gamma _l$, $l=0,\ldots ,q$:

$$\begin{aligned} \dfrac{\partial \ell (\varvec{\theta })}{\partial \beta _j} = \sum _{i=1}^n {\bar{y}}_{ij} \big ( \tau _i - \pi _i \big ), \qquad \dfrac{\partial \ell (\varvec{\theta })}{\partial \gamma _l} = -\sum _{i=1}^n {\bar{w}}_{il} \,\tau _i\, a_i \end{aligned}$$

with $a_i = 1- r_i + (m-1)(1-\xi _i)$. Accordingly, matrix (19) is obtained from the column-by-row product of the incomplete score vector with its transpose.