A nested expectation-maximization algorithm for latent class models with covariates

TL;DR

This paper introduces a nested EM algorithm for latent class models with covariates that ensures monotonic log-likelihood improvement and faster convergence compared to existing methods.

Contribution

It presents a novel nested EM routine that guarantees monotone log-likelihood sequences and enhances convergence speed for latent class models with covariates.

Findings

Guarantees monotone increase in log-likelihood

Achieves improved convergence rates

Applicable to latent class models with covariates

Abstract

We develop a nested EM routine for latent class models with covariates which allows maximization of the full-model log-likelihood and, differently from current methods, guarantees monotone log-likelihood sequences along with improved convergence rates.

Tables

Table 1. Table 1: Quantitative evaluation of the maximization quality and the computational efficiency of the routines discussed in Sect.  1.1 – 2 . The value max 𝜽 ​ { ℓ ​ ( 𝜽 ; 𝐲 , 𝐱 ) } subscript max 𝜽 ℓ 𝜽 𝐲 𝐱 \mbox{max}_{{\boldsymbol{\theta}}}\{\ell({\boldsymbol{\theta}};{\bf y},{\bf x})\} coincides with the highest log-likelihood observed in the 100 runs of the different algorithms. Since most of these routines are devised to maximize ℓ ​ ( 𝜽 ; 𝐲 , 𝐱 ) ℓ 𝜽 𝐲 𝐱 \ell({\boldsymbol{\theta}};{\bf y},{\bf x}) , such choice can be safely regarded as the maximum log-likelihood.

ELECTION DATA $R=3$	nrEM	nrEM$Q_1$	nrEM$Q_1$ $\alpha =0.5$	mmEM
Number of runs with a decay in $\ell ({𝜽}^{\left(t\right)};𝐲,𝐱)$	78	37	12	0
Number of runs reaching a local mode	94	51	24	27
$\left|\ell (\widehat{𝜽};𝐲,𝐱)-{\text{max}}_{𝜽}\left\{\ell (𝜽;𝐲,𝐱)\right\}\right|$ local modes	830.046	1105.626	4.894	0.644
Iterations to reach ${\text{max}}_{𝜽}\left\{\ell (𝜽;𝐲,𝐱)\right\}$	146	151	161	229
Averaged computational time for each run	0.073”	0.182”	0.230”	0.283”
ELECTION DATA $R=3$	3stepClassical	3stepCorrection	nestedEM	hybridEM
Number of runs with a decay in $\ell ({𝜽}^{\left(t\right)};𝐲,𝐱)$	na	na	0	0
Number of runs reaching a local mode	100	100	24	25
$\left|\ell (\widehat{𝜽};𝐲,𝐱)-{\text{max}}_{𝜽}\left\{\ell (𝜽;𝐲,𝐱)\right\}\right|$ local modes	42.224	39.104	0.644	0.644
Iterations to reach ${\text{max}}_{𝜽}\left\{\ell (𝜽;𝐲,𝐱)\right\}$	na	na	171	166
Averaged computational time for each run	0.229”	0.212”	0.359”	0.265”

Equations

$$\text{mbox}pr({\mathbf{Y}}_i=\mathbf{y}\mid {\mathbf{x}}_i)=\sum _{r=1}^R{\nu }_r({\mathbf{x}}_i)\prod _{j=1}^J{\pi }_{jr}(y_j)=\sum _{r=1}^R\frac{\exp ({\mathbf{x}}_i^{⊺}{\mathbf{\beta }}_r)}{{\sum }_{l=1}^R\exp ({\mathbf{x}}_i^{⊺}{\mathbf{\beta }}_l)}\prod _{j=1}^J{\pi }_{jr}(y_j),\text{mbox}foreachi=1,\dots ,n,$$

$$\ell (\mathbf{\beta },\mathbf{\pi };\mathbf{y},\mathbf{x})=\sum _{i=1}^n\log \left[\sum _{r=1}^R{\nu }_r({\mathbf{x}}_i)\prod _{j=1}^J\prod _{y_j=1}^{K_j}{\pi }_{jr}(y_j{)}^{1(y_{ij}=y_j)}\right]=\sum _{i=1}^n\log \left[\sum _{r=1}^R\frac{\exp ({\mathbf{x}}_i^{⊺}{\mathbf{\beta }}_r)}{{\sum }_{l=1}^R\exp ({\mathbf{x}}_i^{⊺}{\mathbf{\beta }}_l)}\prod _{j=1}^J\prod _{y_j=1}^{K_j}{\pi }_{jr}(y_j{)}^{1(y_{ij}=y_j)}\right],$$

$$(Y_{ij}\mid s_i=r)\sim \text{mbox}categorical({\mathbf{\pi }}_{jr},K_j),\text{mbox}foranyj=1,\dots ,J,(s_i\mid {\mathbf{x}}_i)\sim \text{mbox}categorical(\mathbf{\nu }({\mathbf{x}}_i),R),$$

$$\begin{array}{r}{\ell }_1(\mathbf{\beta };\mathbf{s},\mathbf{x})+{\ell }_2(\mathbf{\pi };\mathbf{y},\mathbf{s})=\sum _{i=1}^n\sum _{r=1}^{R}1(s_i=r)\log \left[\frac{\exp ({\mathbf{x}}_i^{⊺}{\mathbf{\beta }}_r)}{{\sum }_{l=1}^R\exp ({\mathbf{x}}_i^{⊺}{\mathbf{\beta }}_l)}\right]+\sum _{i=1}^n\sum _{r=1}^R\left[\sum _{j=1}^J\sum _{y_j=1}^{K_j}1(s_i=r)1(y_{ij}=y_j)\log {\pi }_{jr}(y_j)\right].\end{array}$$

$${\hat{\pi }}_{jr}(y_j)=\frac{{\sum }_{i=1}^{n}1(s_i=r)1(y_{ij}=y_j)}{{\sum }_{i=1}^{n}1(s_i=r),}\text{mbox}foreachr=1,\dots ,R,j=1,\dots ,J,y_j=1,\dots ,K_j.$$

$${\bar{s}}_{ir}^{(t)}=\text{mbox}E\{1(s_i=r)\mid {\mathbf{\theta }}^{(t)},{\mathbf{y}}_i,{\mathbf{x}}_i\}=\frac{\exp \{{\mathbf{x}}_i^{⊺}{\mathbf{\beta }}_r^{(t)}\}{\prod }_{j=1}^J{\prod }_{y_j=1}^{K_j}{\pi }_{jr}^{(t)}(y_j{)}^{1(y_{ij}=y_j)}}{{\sum }_{l=1}^R\exp \{{\mathbf{x}}_i^{⊺}{\mathbf{\beta }}_l^{(t)}\}{\prod }_{j=1}^J{\prod }_{y_j=1}^{K_j}{\pi }_{jl}^{(t)}(y_j{)}^{1(y_{ij}=y_j)}},\text{mbox}foreachr=1,\dots ,R,i=1,\dots ,n,$$

$${\pi }_{jr}^{(t+1)}(y_j)=\frac{{\sum }_{i=1}^n{\bar{s}}_{ir}^{(t)}1(y_{ij}=y_j)}{{\sum }_{i=1}^n{\bar{s}}_{ir}^{(t)}},\text{mbox}foreachr=1,\dots ,R,j=1,\dots ,J,y_j=1,\dots ,K_j.$$

$$Q_1(\mathbf{\beta }\mid {\mathbf{\theta }}^{(t)})=Q_1({\mathbf{\beta }}_1\mid {\mathbf{\theta }}^{(t)})=\sum _{i=1}^n\log \left[\frac{\exp ({\mathbf{x}}_i^{⊺}{\mathbf{\beta }}_1{)}^{{\bar{s}}_{i1}^{(t)}}}{\{1+\exp ({\mathbf{x}}_i^{\mathrm{T}}{\mathbf{\beta }}_1){\}}^{{\bar{s}}_{i1}^{(t)}}}\cdot \frac{\{1+\exp ({\mathbf{x}}_i^{⊺}{\mathbf{\beta }}_1){\}}^{{\bar{s}}_{i1}^{(t)}}}{1+\exp ({\mathbf{x}}_i^{⊺}{\mathbf{\beta }}_1)}\right]=\sum _{i=1}^n\log \left[\frac{\exp ({\mathbf{x}}_i^{⊺}{\mathbf{\beta }}_1{)}^{{\bar{s}}_{i1}^{(t)}}}{1+\exp ({\mathbf{x}}_i^{⊺}{\mathbf{\beta }}_1)}\right],$$

$${\ell }_1^{\ast }({\mathbf{\beta }}_1;{\bar{\mathbf{s}}}^{(t)},\mathbf{x},\mathbf{\omega })$$

$${\mathbf{\beta }}_1^{(t+1)}=({\mathbf{X}}^{⊺}{\bar{\mathbf{\Omega }}}^{(t)}\mathbf{X}{)}^{-1}{\mathbf{X}}^{⊺}{\bar{\mathbf{\Omega }}}^{(t)}{\bar{\mathbf{\eta }}}^{(t)},$$

$$Q_1({\mathbf{\beta }}^{(t+r/R^{\ast })}\mid {\mathbf{\pi }}^{(t+1)},{\mathbf{\beta }}^{\{t+(r-1)/R^{\ast }\}})\ge Q_1({\mathbf{\beta }}^{\{t+(r-1)/R^{\ast }\}}\mid {\mathbf{\pi }}^{(t+1)},{\mathbf{\beta }}^{\{t+(r-1)/R^{\ast }\}}),\text{mbox}foreachr=1,\dots ,R^{\ast }.$$

$$\sum _{i=1}^n\left\{{\bar{s}}_{ir}^{\{t+(r-1)/R^{\ast }\}}\log \left[\frac{\exp ({\mathbf{x}}_i^{⊺}{\mathbf{\beta }}_r)}{\exp ({\mathbf{x}}_i^{⊺}{\mathbf{\beta }}_r)+c_i^{\{t+(r-1)/R^{\ast }\}}}\right]+\sum _{l\ne r}{\bar{s}}_{il}^{\{t+(r-1)/R^{\ast }\}}\log \left[\frac{\exp ({\mathbf{x}}_i^{⊺}{\mathbf{\beta }}_l^{\{t+(r-1)/R^{\ast }\}})}{\exp ({\mathbf{x}}_i^{⊺}{\mathbf{\beta }}_r)+c_i^{\{t+(r-1)/R^{\ast }\}}}\right]\right\},$$

$$Q_1({\mathbf{\beta }}_r\mid {\mathbf{\theta }}^{\{t+(r-1)/R^{\ast }\}})=\sum _{i=1}^n\log \left[\frac{\{\exp ({\mathbf{x}}_i^{⊺}{\mathbf{\beta }}_r-a_i^{\{t+(r-1)/R^{\ast }\}}){\}}^{{\bar{s}}_{ir}^{\{t+(r-1)/R^{\ast }\}}}}{1+\exp ({\mathbf{x}}_i^{⊺}{\mathbf{\beta }}_r-a_i^{\{t+(r-1)/R^{\ast }\}})}\right]+\text{mbox}const,$$

$${\bar{\omega }}_{ir}^{\{t+(r-1)/R^{\ast }\}}=0.5({\mathbf{x}}_i^{⊺}{\mathbf{\beta }}_r^{(t)}-a_i^{\{t+(r-1)/R^{\ast }\}}{)}^{-1}\text{mbox}tanh[0.5({\mathbf{x}}_i^{⊺}{\mathbf{\beta }}_r^{(t)}-a_i^{\{t+(r-1)/R^{\ast }\}})],$$

$${\mathbf{\beta }}_r^{(t+1)}=({\mathbf{X}}^{⊺}{\bar{\mathbf{\Omega }}}^{\{t+(r-1)/R^{\ast }\}}\mathbf{X}{)}^{-1}{\mathbf{X}}^{⊺}{\bar{\mathbf{\Omega }}}^{\{t+(r-1)/R^{\ast }\}}{\bar{\mathbf{\eta }}}^{\{t+(r-1)/R^{\ast }\}},$$

$$\begin{array}{rl}& \ell ({\mathbf{\pi }}^{(t)},{\mathbf{\beta }}_1^{(t)},\dots ,{\mathbf{\beta }}_r^{(t)},\dots ,{\mathbf{\beta }}_{R^{\ast }}^{(t)};\mathbf{y},\mathbf{x})\le \ell ({\mathbf{\pi }}^{(t+1)},{\mathbf{\beta }}_1^{(t)},\dots ,{\mathbf{\beta }}_r^{(t)},\dots ,{\mathbf{\beta }}_{R^{\ast }}^{(t)};\mathbf{y},\mathbf{x})\le \ell ({\mathbf{\pi }}^{(t+1)},{\mathbf{\beta }}_1^{(t+1)},\dots ,{\mathbf{\beta }}_r^{(t)},\dots ,{\mathbf{\beta }}_{R^{\ast }}^{(t)};\mathbf{y},\mathbf{x})\le \\ & \dots \le \ell ({\mathbf{\pi }}^{(t+1)},{\mathbf{\beta }}_1^{(t+1)},\dots ,{\mathbf{\beta }}_r^{(t+1)},\dots ,{\mathbf{\beta }}_{R^{\ast }}^{(t)};\mathbf{y},\mathbf{x})\le \dots \le \ell ({\mathbf{\pi }}^{(t+1)},{\mathbf{\beta }}_1^{(t+1)},\dots ,{\mathbf{\beta }}_r^{(t+1)},\dots ,{\mathbf{\beta }}_{R^{\ast }}^{(t+1)};\mathbf{y},\mathbf{x}).\end{array}$$

Full text

††journal: Statistics and Probability Letters

1 Introduction

Multivariate categorical data are routinely collected in several fields [e.g. 14]. In these settings, it is of key interest to characterize the dependence structures in the observed data, and to identify underlying classes or subpopulations which may explain these patterns of dependence and their changes with external covariates. Let ${\bf Y}_{i}=(Y_{i1},\ldots,Y_{iJ})^{\intercal}$ denote the multivariate categorical random variable generating the observed data ${\bf y}_{i}=(y_{i1},\ldots,y_{iJ})^{\intercal}\in\ \mathcal{Y}=\{1,\ldots,K_{1}\}$ $\times\cdots\times\{1,\ldots,K_{J}\}$, for every unit $i=1,\ldots,n$. Latent class models with covariates [3, 9, 12] address this goal by assuming the response variables $Y_{i1},\ldots,Y_{iJ}$, are conditionally independent given a latent class indicator $s_{i}\in\mathcal{S}=\{1,\ldots,R\}$, whose probability mass function is allowed to change with the covariates ${\bf x}_{i}=(x_{i1},\ldots,x_{iP})^{\intercal}$ $\in\mathcal{X}$, under a multinomial logistic regression. Consistent with this assumption, the conditional probability mass function $\mbox{pr}({\bf Y}_{i}={\bf y}\mid{\bf x}_{i})=\mbox{pr}(Y_{i1}=y_{1},\ldots,Y_{iJ}=y_{J}\mid{\bf x}_{i})$ for the multivariate random variable ${\bf Y}_{i}$, can be expressed as

[TABLE]

where $\nu_{r}({\bf x}_{i})=\mbox{pr}(s_{i}=r\mid{\bf x}_{i})\in(0,1)$ is the covariate–dependent probability of class $r$, whereas $\pi_{jr}(y_{j})=\mbox{pr}(Y_{ij}=y_{j}\mid s_{i}=r)\in(0,1)$ characterizes the probability to observe the category $y_{j}$ for the variable $Y_{ij}$ in class $r$. Note also that, consistent with classical multinomial logistic regression, the coefficients vector $\boldsymbol{\beta}_{R}=(\beta_{1R},\ldots,\beta_{PR})^{\intercal}$ associated with the last class $R$ is fixed to zero, in order to avoid identifiability issues. With this choice, ${\bf x}_{i}^{\intercal}\boldsymbol{\beta}_{r}$ measures the log-odds of belonging to class $r$ instead of class $R$, when the vector of covariates is ${\bf x}_{i}$. Equation (1) provides an interpretable factorization which allows inference on the class-specific generative mechanisms underlying the observed data ${\bf y}_{i}$ and how the latent classes $s_{i}$ relate to the covariates ${\bf x}_{i}$. Refer to [3, 9, 12] for a discussion about this class of models, and to [8] for a review on latent class analysis, including early formulations without covariates [17, 21].

To obtain the above information, it is necessary to estimate $\boldsymbol{\beta}=\{\boldsymbol{\beta}_{1},\ldots,\boldsymbol{\beta}_{R-1}\}$ and $\boldsymbol{\pi}=\{\pi_{j1}(y_{j}),\ldots,\pi_{jR}(y_{j}):j=1,\ldots,J;\ y_{j}=1,\ldots,K_{j}\}$. This can be accomplished by maximizing the log-likelihood function

[TABLE]

where ${1}(y_{ij}=y_{j})$ is $1$ if $y_{ij}=y_{j}$, and [math] otherwise. However, maximization of (2) is not straightforward due to the sum inside the logarithm. In fact, although some contributions attempt direct maximization of (2) via Newton–Raphson [e.g. 13] or simplex algorithms [e.g. 9], more popular implementations [3, 12, 25, 6, 26] rely on em routines [10]. These strategies leverage a hierarchical specification—equivalent to (1)—which introduces a latent class variable $s_{i}\in\mathcal{S}=\{1,\ldots,R\}$ for each $i=1,\ldots,n$, to obtain

[TABLE]

independently for every unit $i=1,\ldots,n$. In (3), $\mbox{{categorical}}(\boldsymbol{\rho},H)$ denotes the generic categorical random variable having probability mass function $\boldsymbol{\rho}=\{{\rho}_{1},\ldots,{\rho}_{H}\}$ for the $H$ different categories. Hence, consistent with (3), if ${\bf s}=(s_{1},\ldots,s_{n})$ is known, the maximum likelihood estimates for $\boldsymbol{\beta}$ and $\boldsymbol{\pi}$ can be easily obtained by maximizing separately the log-likelihood $\ell_{1}(\boldsymbol{\beta};{\bf s},{\bf x})$ associated with the multinomial logistic regression for $s_{1},\ldots,s_{n}$, and the log-likelihood $\ell_{2}(\boldsymbol{\pi};{\bf y},{\bf s})$ for the categorical data $y_{1j},\ldots,y_{nj}$, $j=1,\ldots,J$, within each subpopulation defined by the classes in ${\bf s}$. In fact, $\ell(\boldsymbol{\beta},\boldsymbol{\pi};{\bf y},{\bf s},{\bf x})=\ell_{1}(\boldsymbol{\beta};{\bf s},{\bf x})+\ell_{2}(\boldsymbol{\pi};{\bf y},{\bf s})$, with

[TABLE]

Maximizing $\ell_{1}(\boldsymbol{\beta};{\bf s},{\bf x})$ with respect to $\boldsymbol{\beta}$ requires algorithms for multinomial logistic regression, within a generalized linear models framework [e.g. 1, Chapter 7], whereas $\ell_{2}(\boldsymbol{\pi};{\bf y},{\bf s})$ is analytically maximized at

[TABLE]

Unfortunately, ${\bf s}$ is not observed. Therefore, estimation of $\boldsymbol{\beta}$ and $\boldsymbol{\pi}$ needs to rely only on the information provided by the data $({\bf y}_{i},{\bf x}_{i})$, for $i=1,\ldots,n$.

There are two main strategies in the literature to estimate $\boldsymbol{\beta}$ and $\boldsymbol{\pi}$, generally referred to as one–step [3, 12, 25] and three–step [6, 26] methods. The former attempt simultaneous estimation of $\boldsymbol{\beta}$ and $\boldsymbol{\pi}$ in (2), leveraging the augmented log-likelihood (4). The latter consider, instead, a multi–step routine which first estimates $s_{1},\ldots,s_{n}$, along with $\boldsymbol{\pi}$, from a latent class model without covariates, and then uses the predicted classes as responses in $\ell_{1}(\boldsymbol{\beta};\hat{{\bf s}},{\bf x})$—or a modification of it—to estimate $\boldsymbol{\beta}$.

Although the above methods are considered in routine implementations—including the R library poLCA [19] and the software Latent GOLD [27]—as discussed in Sect. 1.1–1.2, both strategies still raise concerns on the quality of the estimates and on the efficiency of the algorithms. Motivated by these issues, Sect. 2 describes a nested em which provides more reliable estimation routines for this class of models within the maximum likelihood framework. As outlined on a real dataset in Sect. 3, the proposed methods enjoy improved theoretical properties and superior performance. Concluding remarks can be found in Sect. 4.

1.1 One–step estimation methods

Recalling Sect. 1, maximum likelihood estimation for the parameters in the full model (1)—characterizing one–step methods [3, 12, 25]—proceeds via an em algorithm which leverages the complete log-likelihood $\ell(\boldsymbol{\beta},\boldsymbol{\pi};{\bf y},{\bf s},{\bf x})=\ell_{1}(\boldsymbol{\beta};{\bf s},{\bf x})+\ell_{2}(\boldsymbol{\pi};{\bf y},{\bf s})$ (4) for the data $({\bf y}_{i},{\bf x}_{i})$ and the augmented latent class variable $s_{i}$, $i=1,\ldots,n$. This additive structure of $\ell(\boldsymbol{\beta},\boldsymbol{\pi};{\bf y},{\bf s},{\bf x})$ allows separate estimation for the parameters $\boldsymbol{\beta}$ and $\boldsymbol{\pi}$. Moreover, $\ell_{1}(\boldsymbol{\beta};{\bf s},{\bf x})$ and $\ell_{2}(\boldsymbol{\pi};{\bf y},{\bf s})$ are linear in the augmented data ${1}(s_{i}=r)$. Letting $\boldsymbol{\theta}=(\boldsymbol{\beta},\boldsymbol{\pi})$, this facilitates a simple expectation step in which, at the general iteration $t$, each ${1}(s_{i}=r)$ is replaced with

[TABLE]

to obtain the expected values $Q_{1}(\boldsymbol{\beta}\mid\boldsymbol{\theta}^{(t)})$ and $Q_{2}(\boldsymbol{\pi}\mid\boldsymbol{\theta}^{(t)})$ of the log-likelihoods $\ell_{1}(\boldsymbol{\beta};{\bf s},{\bf x})$ and $\ell_{2}(\boldsymbol{\pi};{\bf y},{\bf s})$, respectively, whose summation defines the expectation $Q(\boldsymbol{\beta},\boldsymbol{\pi}\mid\boldsymbol{\theta}^{(t)})$ of the complete log-likelihood (4) with respect to the distribution of ${\bf s}$, given the current estimates $\boldsymbol{\theta}^{(t)}$. Hence, $\boldsymbol{\theta}^{(t+1)}=\mbox{argmax}_{\boldsymbol{\beta},\boldsymbol{\pi}}\{Q(\boldsymbol{\beta},\boldsymbol{\pi}\mid\boldsymbol{\theta}^{(t)})\}$ in the m–step can be obtained by maximizing $Q_{1}(\boldsymbol{\beta}\mid\boldsymbol{\theta}^{(t)})$ and $Q_{2}(\boldsymbol{\pi}\mid\boldsymbol{\theta}^{(t)})$ separately.

Consistent with (5), the expected log-likelihood $Q_{2}(\boldsymbol{\pi}\mid\boldsymbol{\theta}^{(t)})$ is easily maximized at

[TABLE]

It is instead not possibile to maximize analytically $Q_{1}(\boldsymbol{\beta}\mid\boldsymbol{\theta}^{(t)})$ with respect to $\boldsymbol{\beta}$, due to the logistic link. To address this issue [12] and [25] consider one Newton–Raphson step relying on a quadratic approximation of $Q_{1}(\boldsymbol{\beta}\mid\boldsymbol{\theta}^{(t)})$. The routine proposed by [3] leverages instead a slightly different update relying on the Hessian and the gradient computed from full–model log-likelihood in (2), and evaluated at $(\boldsymbol{\beta}^{(t)},\boldsymbol{\pi}^{(t)})$ and $\bar{s}^{(t)}_{ir}$, $r=1,\ldots,R$, $i=1,\ldots,n$. This procedure—implemented in R library poLCA [19]—breaks the em rationale, since the m–step for $\boldsymbol{\beta}$ maximizes a plug–in estimate of a quadratic approximation for $\ell(\boldsymbol{\beta};{\bf y},{\bf x})$, instead of $Q_{1}(\boldsymbol{\beta}\mid\boldsymbol{\theta}^{(t)})$. This may lead to less stable behaviors.

Although the above solutions provide standard approaches to obtain $\boldsymbol{\beta}^{(t+1)}$, the resulting computational routines guarantee neither $Q_{1}(\boldsymbol{\beta}^{(t+1)}\mid\boldsymbol{\theta}^{(t)})\geq Q_{1}(\boldsymbol{\beta}\mid\boldsymbol{\theta}^{(t)})$, nor $Q_{1}(\boldsymbol{\beta}^{(t+1)}\mid\boldsymbol{\theta}^{(t)})\geq Q_{1}(\boldsymbol{\beta}^{(t)}\mid\boldsymbol{\theta}^{(t)})$ [4, 5], and therefore the proposed methods are neither an em, nor a generalized em algorithm, respectively [10]. Failure to improve the expected log-likelihood may also affect the monotonicity of the sequence $\ell(\boldsymbol{\beta}^{(t)},\boldsymbol{\pi}^{(t)};{\bf y},{\bf x})$, thereby providing routines which do not meet the basic properties of em, and may not guarantee reliable convergence [e.g. 22, Chapter 1.5.5]. As we will outline in Sect. 3, this issue is not just found in pathological scenarios, but arises also in routine applications, and may substantially affect convergence to the maximum log-likelihood. Although careful routines can be designed to overcome this issue, we shall emphasize that $Q_{1}(\boldsymbol{\beta}\mid\boldsymbol{\theta}^{(t)})$ is defined on a set of latent responses whose expectation (6) changes at every iteration of the algorithm. This setting is more problematic than classical multinomial logit with observed response data, and is further complicated by the need to estimate $\boldsymbol{\pi}$ along with $\boldsymbol{\beta}$ in (1). Hence, unstable updating steps may have major effects on estimation. Even Monte Carlo em routines [28] do not address the problem, since such methods are devised for intractable e–steps and not for m–steps having no analytical solutions.

To mitigate the above issues, standard implementations consider multiple runs based on different initializations of the em, and rely on the routine converging to the highest log-likelihood. Alternatively, internal checks can be included to control for decays. These strategies provide more reliable routines, but require multiple runs which increase computational costs. A possibility to reduce decays in $\ell(\boldsymbol{\beta}^{(t)},\boldsymbol{\pi}^{(t)};{\bf y},{\bf x})$, without relying on multiple runs, is to rescale the inverse of the Hessian by a step–size $0<\alpha\leq 1$ [e.g. 22, Chapter 1.5.6]. However, the implementation requires the choice of $\alpha$, without theory on optimal settings.

Motivated by the above issues, [4, 5] proposed a Minorize–Majorize (mm) algorithm [e.g. 16] for logistic and multinomial logit regression, which replaces the Hessian with a matrix $\boldsymbol{B}$ to obtain a quadratic function minorizing the log-likelihood at every $\boldsymbol{\beta}$ and tangent to it in $\boldsymbol{\beta}^{(t)}$. This provides a simple updating scheme similar to the Newton–Raphson which guarantees monotone log-likelihood sequences. Although not currently implemented in latent class models, these strategies can be easily incorporated in the m–step for $\boldsymbol{\beta}$ by replacing the observed responses with the expectation (6) of the latent classes. As outlined in Sect.3, this procedure provides monotone log-likelihood sequences, but requires more iterations than nested em.

1.2 Three–step estimation methods

Differently from one–step methods, the three–step procedures [e.g. 8] do not attempt direct maximization of (2), but instead rely on a multi–step strategy which first estimates $\boldsymbol{\pi}$ and the class probabilities via closed–form em for latent class models without covariates, and then obtain $\hat{\boldsymbol{\beta}}$ from a multinomial logit with the predicted classes $\hat{s}_{1},\ldots,\hat{s}_{n}$ from the previous model acting as observed responses.

The above procedure provides a simple estimation strategy, which is motivated by its interpretability. However, as discussed in [6] and [26], the resulting estimators for the coefficients in $\boldsymbol{\beta}$ are subject to systematic bias. Motivated by these issues, [6] and [26], recently developed two bias–correction methods relying on a modification of the multinomial log-likelihood for $\hat{s}_{1},\ldots,\hat{s}_{n}$, which incorporates the classification error. The bias–corrected log-likelihood proposed in [6] allows simple estimation via standard algorithms for multinomial logistic regression, but this method can be applied only when the covariates are categorical. A more general correction procedure is considered by [26] to include continuous covariates. This approach further improves efficiency and facilitates wider applicability. There is also a focus on other mechanisms to predict the classes, but the final results are not substantially different [26, 2].

Although the bias adjustment in [26] is widely considered in several implementations, the estimates remain still sub-optimal compared to one–step methods, since they do not directly maximize the full–model log-likelihood (2), and hence cannot be considered as the maximum likelihood estimates. Indeed, when the focus is on providing reliable inference for the parameters $\boldsymbol{\beta}$ and $\boldsymbol{\pi}$ in (1), it is arguably more coherent to attempt a direct maximization of the log-likelihood in (2), since it guarantees unbiased, efficient and consistent estimators, under the classical likelihood inference theory, if the model is correctly specified [6].

2 Nested em for one–step estimation

To address the aforementioned issues, we propose a nested em algorithm for one–step estimation which avoids approximations of the expected log-likelihood $Q_{1}(\boldsymbol{\beta}\mid\boldsymbol{\theta}^{(t)})$, but improves this function sequentially via a set of conditional expectation–maximizations for every vector of coefficients $\boldsymbol{\beta}_{r}$, given the others.

Working with conditional expected log-likelihoods is appealing in providing a set of different logistic regressions for which the recent Pólya-gamma data augmentation scheme [24, 7] guarantees closed–form maximization via generalized least squares. Indeed, following Theorem 1 in [24], the generic logistic likelihood $\exp({\bf x}^{\intercal}\boldsymbol{\beta})^{a}\{1+\exp({\bf x}^{\intercal}\boldsymbol{\beta})\}^{-b}$, can be rewritten as $2^{-b}\exp\{(a-0.5b){\bf x}^{\intercal}\boldsymbol{\beta}\}\mbox{cosh}(0.5{\bf x}^{\intercal}\boldsymbol{\beta})^{-b}$, $b\geq 0$, whereas the likelihood of a Pólya-gamma variable $\omega\sim\textsc{pg}(b,{\bf x}^{\intercal}\boldsymbol{\beta})$ is proportional to $\exp\{-0.5\omega({\bf x}^{\intercal}\boldsymbol{\beta})^{2}\}\mbox{cosh}(0.5{\bf x}^{\intercal}\boldsymbol{\beta})^{b}$. Hence, combining these two quantities provides an augmented likelihood $\exp\{-0.5\omega({\bf x}^{\intercal}\boldsymbol{\beta})^{2}+(a-0.5b){\bf x}^{\intercal}\boldsymbol{\beta}\}$ for data $(a,\omega)$ which is proportional to the one induced by a Gaussian regression for the transformed response $\omega^{-1}(a-0.5b)\sim\mbox{N}({\bf x}^{\intercal}\boldsymbol{\beta},\omega^{-1})$. This allows simple maximization for $\boldsymbol{\beta}$ within a Gaussian framework. Note that, as discussed in [24], other Gaussian–related data augmentations for logistic regression have been developed. However, these strategies require more complex representations [24], thus leading to intractable computations. Instead, the Pólya-gamma data augmentation leads directly to a Gaussian likelihood with a single latent variable $\omega$, whose expectation is analytically available via $\mbox{E}(\omega)=0.5b({\bf x}^{\intercal}\boldsymbol{\beta})^{-1}\mbox{tanh}(0.5{\bf x}^{\intercal}\boldsymbol{\beta})$. Besides providing tractable computations, the proposed procedure guarantees monotone log-likelihood sequences and is directly motivated by an exact em for the special case of $R=2$ classes, described below.

2.1 Exact em algorithm for $R=2$

Let us first focus on deriving an em for $R=2$ latent classes, which provides analytical maximization also for $\boldsymbol{\beta}_{1}$. In fact, when $R=2$ the expected value of $\ell_{1}(\boldsymbol{\beta};{\bf s},{\bf x})=\ell_{1}(\boldsymbol{\beta}_{1};{\bf s},{\bf x})$ is

[TABLE]

thus providing the log-likelihood $\ell^{*}_{1}(\boldsymbol{\beta}_{1};\bar{\bf s}^{(t)},{\bf x})$ of a logistic regression in which the entries in $\bar{\bf s}^{(t)}$ act as responses. This allows the implementation of the Pólya-gamma data augmentation. In particular, defining $b:=1$, $a:=\bar{s}_{i1}^{(t)}$, ${\bf x}^{\intercal}\boldsymbol{\beta}:={\bf x}_{i}^{\intercal}\boldsymbol{\beta}_{1}$ and $\omega:=\omega_{i1}$, leads to the complete log-likelihood

[TABLE]

for $\bar{\bf s}^{(t)}$, ${\bf x}$, and the Pólya-gamma augmented data $\boldsymbol{\omega}=(\omega_{11},\ldots,\omega_{n1})$. Equation (8) is a quadratic function of ${\bf x}_{i}^{\intercal}\boldsymbol{\beta}_{1}$, and is linear in $\boldsymbol{\omega}$. This allows the implementation of a simple nested expectation step in which every $\omega_{i1}$ is replaced with the expectation $\mbox{E}(\omega_{i1}\mid\boldsymbol{\pi}^{(t)},\boldsymbol{\beta}_{1}^{(t)},{\bf y}_{i},{\bf x}_{i}):=\bar{\omega}^{(t)}_{i1}=0.5({\bf x}_{i}^{\intercal}\boldsymbol{\beta}^{(t)}_{1})^{-1}\mbox{tanh}(0.5{\bf x}_{i}^{\intercal}\boldsymbol{\beta}^{(t)}_{1})$ to obtain $Q^{*}_{1}(\boldsymbol{\beta}_{1}\mid\boldsymbol{\theta}^{(t)})=\sum_{i=1}^{n}-0.5\bar{\omega}^{(t)}_{i1}(\bar{\eta}^{(t)}_{i1}-{\bf x}_{i}^{\intercal}\boldsymbol{\beta}_{1})^{2}+\mbox{const}$, with $\bar{\eta}^{(t)}_{i1}=(\bar{s}_{i1}^{(t)}-0.5)/\bar{\omega}^{(t)}_{i1}$. The appealing property associated with this nested expected log-likelihood, compared to $Q_{1}(\boldsymbol{\beta}_{1}\mid\boldsymbol{\theta}^{(t)})$, is that it allows direct maximization for $\boldsymbol{\beta}_{1}$. Indeed, exploiting the generalized least squares, $Q^{*}_{1}(\boldsymbol{\beta}_{1}\mid\boldsymbol{\theta}^{(t)})$ is maximized at

[TABLE]

where ${\bf X}$ is the $n\times P$ matrix having rows ${\bf x}^{\intercal}_{i}$, whereas $\bar{\boldsymbol{\Omega}}^{(t)}=\mbox{diag}(\bar{\omega}^{(t)}_{11},\ldots,\bar{\omega}^{(t)}_{n1})$ and $\bar{\boldsymbol{\eta}}^{(t)}=(\bar{\eta}^{(t)}_{11},\ldots,\bar{\eta}^{(t)}_{n1})^{\intercal}$. Refer to [11] for recent results proving that solution (9) guarantees monotone convergence at optimal rate in logistic regression. This provides additional support for (9), compared to Newton–Raphson and mm updates.

We shall stress that, although we obtain $Q^{*}_{1}(\boldsymbol{\beta}_{1}\mid\boldsymbol{\theta}^{(t)})$ sequentially, this quantity coincides with the expectation of the complete log-likelihood $\ell_{1}(\boldsymbol{\beta}_{1};{\bf s},{\bf x},\boldsymbol{\omega})$. Hence, since $Q_{2}(\boldsymbol{\pi}\mid\boldsymbol{\theta}^{(t)})$ is analytically maximized in (7), the resulting routine is an exact em algorithm based on the complete log-likelihood $\ell(\boldsymbol{\beta}_{1},\boldsymbol{\pi};{\bf y},{\bf s},{\bf x},\boldsymbol{\omega})=\ell_{1}(\boldsymbol{\beta}_{1};{\bf s},{\bf x},\boldsymbol{\omega})+\ell_{2}(\boldsymbol{\pi};{\bf y},{\bf s})$. Finally, note that since $\boldsymbol{\pi}^{(t+1)}$ can be easily obtained before updating $\boldsymbol{\beta}$, a more efficient strategy, inspired by the multi-cycle expectation conditional maximization algorithm [23], is to update also the expectation of ${\bf s}$ based on $\boldsymbol{\pi}^{(t+1)}$ instead of $\boldsymbol{\pi}^{(t)}$, before applying (9). As discussed in Sect. 2.2, this solution will be adopted in the general case with $R>2$.

2.2 Nested em algorithm for $R>2$

When more than two classes are considered, $Q_{1}(\boldsymbol{\beta}\mid\boldsymbol{\theta}^{(t)})$ has a multinomial logit form, and not a logistic one. Thus, a direct application of the Pólya-gamma data augmentation is not possible. However, as we will outline, the conditional expected log-likelihood for every vector of coefficients $\boldsymbol{\beta}_{r}$, given the others, can be rewritten as a proper logistic log-likelihood [e.g. 15], thus motivating a Pólya-gamma data augmentation. Based on this result, we propose a nested em which improves the expected log-likelihood for $\boldsymbol{\beta}$ via a set of conditional expectation–maximizations. In particular, for each iteration $t$, we consider $R^{*}=R-1$ nested cycles which sequentially improve the conditional expected log-likelihood of one vector of coefficients $\boldsymbol{\beta}_{r}$, fixing the other $\boldsymbol{\beta}_{l}$, $l\neq r$ at their most recent value. Hence, let $\boldsymbol{\beta}^{(t+r/R^{*})}=\{\boldsymbol{\beta}_{1}^{(t+1)},\ldots,\boldsymbol{\beta}^{(t+1)}_{r},\boldsymbol{\beta}^{(t)}_{r+1},\ldots,\boldsymbol{\beta}^{(t)}_{R^{*}}\},$ denote the estimates for the class-specific vectors of coefficients at cycle $r=1,\ldots,R^{*}$ in iteration $t$, we seek a sequential updating procedure providing the chain inequalities

[TABLE]

The key difference between $\boldsymbol{\beta}^{(t+r/R^{*})}$ and $\boldsymbol{\beta}^{\{t+(r-1)/R^{*}\}}$ in (10), is that only the coefficients in $\boldsymbol{\beta}_{r}$ are updated, whereas all the others are kept fixed. Hence, at every cycle $r$ we seek to improve the expected log-likelihood produced in $r-1$, by modifying only $\boldsymbol{\beta}_{r}$ from its previous estimate at $t$ to a new one at $t+1$. In this respect, such strategy partially recalls the univariate version of Newton–Raphson in [27], and adapts it to blocks of parameters to obtain more efficient updates relying on the most recent estimates. Moreover, as outlined in Proposition 2.1, the sequential improvements underlying the nested em, along with the direct maximization of $Q_{2}(\boldsymbol{\pi}\mid\boldsymbol{\theta}^{(t)})$, guarantee the monotonicity of the sequence $\ell(\boldsymbol{\beta}^{(t)},\boldsymbol{\pi}^{(t)};{\bf y},{\bf x})$. This is fundamental to ensure reliable convergence [4]. Note that in (10), also the expectation of the conditional log-likelihood with respect to the augmented data is sequentially updated using the estimates of the coefficients from the previous cycle, in the same spirit of the multi-cycle expectation conditional maximization [23].

In deriving an updating procedure for $\boldsymbol{\beta}$ with property (10), we adapt results in Sect. 2.1, which provides simple and explicit maximization for each $\boldsymbol{\beta}_{r}$. Focusing on cycle $r$ within iteration $t$, let $Q_{1}(\boldsymbol{\beta}_{r}\mid\boldsymbol{\theta}^{\{t+(r-1)/R^{*}\}})$, with $\boldsymbol{\theta}^{\{t+(r-1)/R^{*}\}}=(\boldsymbol{\pi}^{(t+1)},\boldsymbol{\beta}^{\{t+(r-1)/R^{*}\}})$, denote the conditional expected log-likelihood, written as a function of only $\boldsymbol{\beta}_{r}$, with all the other class-specific coefficients fixed at their corresponding estimates at cycle $r-1$. According to equations (4) and (6), the function $Q_{1}(\boldsymbol{\beta}_{r}\mid\boldsymbol{\theta}^{\{t+(r-1)/R^{*}\}})$ can be expressed as

[TABLE]

where the constants $c_{i}^{\{t+(r-1)/R^{*}\}}$ denote the sum of all the exponential quantities $\exp({\bf x}_{i}^{\intercal}\boldsymbol{\beta}_{l})$, $l\neq r$, written as a function of the current estimates for the class-specific coefficients at cycle $r-1$, whereas the expectation of the augmented latent class indicators are calculated in (6) as a function of $\boldsymbol{\pi}^{(t+1)}$ and $\boldsymbol{\beta}^{\{t+(r-1)/R^{*}\}}$.

Since we aim to improve $Q_{1}(\boldsymbol{\beta}^{(t)}_{r}\mid\boldsymbol{\theta}^{\{t+(r-1)/R^{*}\}})=Q_{1}(\boldsymbol{\beta}^{\{t+(r-1)/R^{*}\}}\mid\boldsymbol{\theta}^{\{t+(r-1)/R^{*}\}})$ by updating the current estimate of $\boldsymbol{\beta}_{r}$ under the Pólya-gamma data augmentation outlined in the previous section, let us highlight a logistic log-likelihood in (11). Indeed, holding out additive constants not depending on $\boldsymbol{\beta}_{r}$, and dividing both the numerator and the denominator of the arguments in the logarithmic functions by the quantities in the vector ${\bf c}^{\{t+(r-1)/R^{*}\}}$, we easily obtain

[TABLE]

where $a_{i}^{\{t+(r-1)/R^{*}\}}=\log c_{i}^{\{t+(r-1)/R^{*}\}},$ and provided that $\bar{s}^{\{t+(r-1)/\bar{R}\}}_{ir}+\sum_{l\neq r}\bar{s}^{\{t+(r-1)/\bar{R}\}}_{il}=1$. Hence, up to the constants $a_{i}^{\{t+(r-1)/\bar{R}\}}$ in the linear predictor, equation (12) for $\boldsymbol{\beta}_{r}$ has the same form of the expected log-likelihood for $\boldsymbol{\beta}_{1}$ in Sect. 2.1, thereby motivating the Pólya-gamma data augmentation at each cycle $r$ of the nested em. In particular, introducing $\omega_{ir}\ {\sim}\ \textsc{pg}(1,{\bf x}_{i}^{\intercal}\boldsymbol{\beta}_{r}-a_{i}^{\{t+(r-1)/R^{*}\}})$, for every $i=1,\ldots,n$, we can exploit the analytical results in Sect. 2.1, to show that the conditional expectation of each $\omega_{ir}$ is

[TABLE]

and that the desired increment (10) at cycle $r$, can be simply obtained—similarly to (9)—by setting

[TABLE]

where $\bar{\boldsymbol{\Omega}}^{\{t+(r{-}1)/R^{*}\}}$ is an $n\times n$ diagonal matrix with elements $\bar{\omega}^{\{t+(r-1)/R^{*}\}}_{ir}$, whereas $\bar{\boldsymbol{\eta}}^{\{t+(r-1)/R^{*}\}}$ is a vector of length $n$ with entries $\bar{\eta}^{\{t+(r-1)/R^{*}\}}_{ir}=(\bar{s}^{\{t+(r-1)/R^{*}\}}_{ir}-0.5+\bar{\omega}^{\{t+(r-1)/R^{*}\}}_{ir}a^{\{t+(r-1)/R^{*}\}}_{i})/\bar{\omega}^{\{t+(r-1)/R^{*}\}}_{ir}$, for $i=1,\ldots,n$.

Algorithm 1 provides details for the implementation of nested em. Note that all the steps require simple and exact expressions, in contrast with current routines. Moreover, as proved in Proposition 2.1, our methods imply monotone convergence.

Proposition 2.1**.**

The nested em in Sect. 2.2 implies $\ell(\boldsymbol{\beta}^{(t+1)},\boldsymbol{\pi}^{(t+1)};{\bf y},{\bf x})\geq\ell(\boldsymbol{\beta}^{(t)},\boldsymbol{\pi}^{(t)};{\bf y},{\bf x})$ at any $t$.

Proof.

To prove Proposition 2.1, we start by showing that the nested conditional expectation–maximizations for each $\boldsymbol{\beta}_{r}$, $r=1,\ldots,R^{*}$, imply (10). Since each of these optimizations steps is a pure em based on the Pólya-gamma data augmentation, following [10], we have $\ell^{*}_{1}(\boldsymbol{\beta}^{(t+1)}_{r};\bar{\bf s}^{\{t+(r-1)/R^{*}\}}_{r},{\bf x})\geq\ell^{*}_{1}(\boldsymbol{\beta}^{(t)}_{r};\bar{\bf s}^{\{t+(r-1)/R^{*}\}}_{r},{\bf x})$ for each $r=1,\ldots,R^{*}$. This implies $Q_{1}(\boldsymbol{\beta}^{\{t+r/R^{*}\}}{\mid}\boldsymbol{\theta}^{\{t+(r-1)/R^{*}\}})\geq Q_{1}(\boldsymbol{\beta}^{\{t+(r-1)/R^{*}\}}\mid\boldsymbol{\theta}^{\{t+(r-1)/R^{*}\}})$, provided that in (12) these log-likelihoods differ by an additive constant which does not depend on $\boldsymbol{\beta}_{r}$. To conclude the proof, we need to ensure that the inequalities in (10), along with the direct maximization of $Q_{2}(\boldsymbol{\pi}\mid\boldsymbol{\theta}^{(t)})$, guarantee $\ell(\boldsymbol{\beta}^{(t+1)},\boldsymbol{\pi}^{(t+1)};{\bf y},{\bf x})\geq\ell(\boldsymbol{\beta}^{(t)},\boldsymbol{\pi}^{(t)};{\bf y},{\bf x})$. Let $Q(\boldsymbol{\beta},\boldsymbol{\pi}\mid\boldsymbol{\theta}^{(t)})=Q_{1}(\boldsymbol{\beta}\mid\boldsymbol{\theta}^{(t)})+Q_{2}(\boldsymbol{\pi}\mid\boldsymbol{\theta}^{(t)})$ be the expected log-likelihood, written as a function of all the parameters in our model. Direct maximization of $Q_{2}(\boldsymbol{\pi}\mid\boldsymbol{\theta}^{(t)})$ in the first step of our routine, guarantees $Q(\boldsymbol{\pi}^{(t)},\boldsymbol{\beta}^{(t)}\mid\boldsymbol{\pi}^{(t)},\boldsymbol{\beta}^{(t)})\leq Q(\boldsymbol{\pi}^{(t+1)},\boldsymbol{\beta}^{(t)}\mid\boldsymbol{\pi}^{(t)},\boldsymbol{\beta}^{(t)})$. Therefore, as discussed in page 165 of [22], this first result guarantees $\ell(\boldsymbol{\pi}^{(t+1)},\boldsymbol{\beta}^{(t)};{\bf y},{\bf x})\geq\ell(\boldsymbol{\pi}^{(t)},\boldsymbol{\beta}^{(t)};{\bf y},{\bf x})$. Similarly, the inequalities in (10), characterizing each nested cycle $r$, ensure $\ell(\boldsymbol{\pi}^{(t+1)},\boldsymbol{\beta}^{(t+r/R^{*})};{\bf y},{\bf x})\geq\ell(\boldsymbol{\pi}^{(t+1)},\boldsymbol{\beta}^{(t+(r-1)/R^{*})};{\bf y},{\bf x})$ for every $r=1,\ldots,R^{*}$. Joining these results we can prove Proposition 2.1 via

[TABLE]

∎

2.3 Hybrid nested em for one–step estimation

Before evaluating the empirical performance of nested em, we first propose a more practical and efficient hybrid modification. Indeed, as discussed in [26], the em is characterized by stable maximization even when the initialization is far from the optimal solution, whereas Newton–Raphson guarantees fast convergence when the routine is close to the maximum. Motivated by this result, we propose an hybrid procedure which starts with nested em and then switches to classical Newton–Raphson [12, 25] when the increment $\ell(\boldsymbol{\beta}^{(t+1)},\boldsymbol{\pi}^{(t+1)};{\bf y},{\bf x})-\ell(\boldsymbol{\beta}^{(t)},\boldsymbol{\pi}^{(t)};{\bf y},{\bf x})$ is less than or equal to a pre–specified small $\epsilon\geq 0$. Although the inclusion of Newton–Raphson steps could still cause decays in the log-likelihood sequence, when the routine is in a neighborhood of the global—or local—maximum, the parabolic approximation is more stable.

3 Empirical study

To study the benefits of nested em in a real–data application, we compare its computational performance with those of the routines discussed in Sect. 1.1 and 1.2. In particular, one–step competitors comprise the em with one Newton–Raphson step (nrEM) from [3], the more formal em (nrEM$Q_{1}$) maximizing $Q_{1}(\boldsymbol{\beta}\mid\boldsymbol{\theta}^{(t)})$ as in [12, 25], and the conservative version of nrEM$Q_{1}$ which relies on the rescaled updating of $\boldsymbol{\beta}$, as outlined in Sect. 1.1, with $\alpha=0.5$ (nrEM$Q_{1}$, $\alpha=0.5$). Recalling Sect. 1.1, we additionally adapt the mm (mmEM) proposed by [4, 5] to compare empirical performance with a routine guaranteeing monotone log-likelihood sequences. The three–step methods considered are instead the classical strategy (3stepClassical) discussed in Sect. 1.2 [e.g. 8], and the popular bias–corrected algorithm (3stepCorrection) from [26]. The hybridEM is instead the modification of our nestedEM presented in Sect. 2.3. In performing maximization under this modified version of the nested em we set $\epsilon=0.01$. Although we have found the results robust to moderate changes in small values of $\epsilon$, it is important to notice that high $\epsilon$ should be avoided since it allows the hybrid routine to rely on Newton–Raphson steps even when $\ell(\boldsymbol{\beta}^{(t)},\boldsymbol{\pi}^{(t)};{\bf y},{\bf x})$ is far from the maximum.

To provide a detailed assessment, we perform estimation under the above algorithms for 100 runs, with varying random initialization. The routines are all initialized at the same values—for each run—and stop when the increment in the log-likelihood is lower than $10^{-11}$. For every run we study the maximization performance and the computational efficiency of the different algorithms. Specifically, the maximization performance is monitored by the number of runs with a decay in the log-likelihood sequence, and by the frequency of runs converging to local modes. For the runs reaching local modes, we also compute the median of the absolute difference between the log-likelihood in these local modes and the maximum one. Computational efficiency is instead studied via the median number of iterations for convergence, computed only for the runs reaching ${\mbox{max}_{{\boldsymbol{\beta,\pi}}}}\{\ell({\boldsymbol{\beta,\pi}};{\bf y},{\bf x})\}$, and the averaged computational time.111Computations rely on R (version 3.3.2) implementations in a machine with 1 Intel Core i5 2.5 GHz processor and 4 GB ram. Code to reproduce the analyses and further results are available at https://github.com/danieledurante/nEM.

Table 1 summarizes the performance of the methods discussed in Sect. 1.1–2, for an application to the election data available in the R library poLCA [19]. This dataset measures voters political affiliation, along with their opinions on how well six different personality traits describe the candidates Al Gore and George Bush before the 2000 presidential elections. These $J=12$ categorical opinions are collected on a four items scale for $n=880$ voters. Here, we assess performance of the maximization routines considering $R=3$ classes, with the political affiliation entering as covariate in the multinomial logit for such latent classes.

Consistent with Sect. 1.1–2, nestedEM and mmEM always provide monotone sequences for (2), thereby guaranteeing accurate maximization and reduced frequency of local modes. Including a Newton–Raphson step as in [3, 12, 25], leads instead to decays in the log-likelihood sequences, which increase the chance of local modes far from ${\mbox{max}_{{\boldsymbol{\beta},\boldsymbol{\pi}}}}\{\ell({\boldsymbol{\beta},\boldsymbol{\pi}};{\bf y},{\bf x})\}$. As discussed in Sect. 1.1, this issue is more severe for nrEM [3] compared to the formal nrEM$Q_{1}$ [12, 25]. This reduced maximization performance of the Newton–type algorithms is mitigated by partial improvements in computational efficiency compared to nestedEM and mmEM, which remain, however, on a similar scale both for the averaged computational time of each run, and for the number of iterations to reach ${\mbox{max}_{{\boldsymbol{\beta},\boldsymbol{\pi}}}}\{\ell({\boldsymbol{\beta},\boldsymbol{\pi}};{\bf y},{\bf x})\}$. Since these algorithms perform estimation in fractions of seconds, an improved maximization performance is arguably the most important property.

Combining nestedEM and nrEM$Q_{1}$, provides an hybridEM which guarantees the accurate maximization performance of monotone algorithms and a computational efficiency comparable to Newton–type methods. Also rescaling the Newton–Raphson updating of $\boldsymbol{\beta}$ by $\alpha=0.5$, allows improvements in maximization performance, but these gains are associated with a reduced computational efficiency. To conclude the evaluation of the one–step algorithms, we shall notice that the correction proposed by [4, 5] guarantees reliable maximization, but, as theoretically proved in [11], their global conservative bound requires more iterations to reach convergence compared to nestedEM and hybridEM.

As discussed in Sect. 1.2, the three–step algorithms do not attempt direct maximization of (2). The consequences of this are evident in Table 1 with all the runs of the 3stepClassical and the 3stepCorrection routines providing sub-optimal estimates which induce local modes in the full–model log-likelihood. Due to this, it is not possible to study the number of iterations to reach ${\mbox{max}_{{\boldsymbol{\beta},\boldsymbol{\pi}}}}\{\ell({\boldsymbol{\beta},\boldsymbol{\pi}};{\bf y},{\bf x})\}$. Also the number of decays in $\ell({\boldsymbol{\beta}}^{(t)},{\boldsymbol{\pi}}^{(t)};{\bf y},{\bf x})$ is somewhat irrelevant to evaluate the three–step methods, since the estimation routines are based on two separate maximizations not directly related to (2). It is also worth noticing that, although our studies are based on the log-likelihood sequence instead of the parameters estimates, within a maximum likelihood approach these two quantities are directly related. In fact, our goal is not on providing a new class of models and estimators for the parameters in (1), but an improved maximization routine for likelihood–based inference, with more reliable convergence to the optimal parameters estimates maximizing (2). In this respect, we have also inspected the estimates associated with the local modes in Table 1 observing evident deviations from the optimal ones, especially in local maxima with a log-likelihood notably below ${\mbox{max}_{{\boldsymbol{\beta},\boldsymbol{\pi}}}}\{\ell({\boldsymbol{\beta},\boldsymbol{\pi}};{\bf y},{\bf x})\}$. This result provides further support to our nested em.

In performing the above studies we initialized $\boldsymbol{\beta}$ from independent Gaussians with mean [math] and small variance $0.5$. Reducing such variance—i.e. initializing $\boldsymbol{\beta}$ close to [math]—yield to gains in all routines. Also in these cases, however, the performance of nestedEM and hybridEM was not worse than the other algorithms, while being less sensitive to initialization. Similar conclusions were obtained in other applications.

4 Discussion

Motivated by the recent Pólya-gamma data augmentation for logistic regression, and by the lack of fully reliable estimation routines for latent class models with covariates, we developed a nested em for one–step maximum likelihood estimation within this class of models. Differently from Newton–type methods, the nested em has theoretical guarantee of monotone log-likelihood sequences, while ensuring convergence at a faster rate than mm routines [11]. This provides improved estimation and computational efficiency.

Although our focus has been on latent class analysis with covariates, the nested em can be also adapted to multinomial logit models with Gaussian random effects. In this setting, the calculation of the expected log-likelihood for the fixed coefficients requires intractable marginalization over the Gaussian random effects, thus requiring Monte Carlo em [28], h–likelihood [18] and other methods [e.g. 20]. Our nested em could provide key benefits over these alternatives. In fact, conditioned on Pólya-gamma augmented data, it is possible to perform closed–form marginalization and maximization as in Gaussian linear mixed models.