The Orthogonally Partitioned EM Algorithm: Extending the EM Algorithm for Algorithmic Stability and Bias Correction Due to Imperfect Data

2.1 The EM algorithm

When there are imperfect data (i.e. missing and mismeasured data), we can obtain MLEs by maximizing the log-likelihood of the observed portion of the data. The log-likelihood associated with the observed data, ℓo, is

where Yi denotes the vector of variables intended for collection from the ith observational unit, f(Yi|θY) is the pdf of Yi parameterized by θY, and the variable Yi={Yi,o,Yi,m} such that Yi,o refers to the elements of Yi that are observable and Yi,m refers to the elements of Yi that are missing. We let yi={yi,o,yi,m} denote the realization of Yi={Yi,o,Yi,m} where yi,o denote the subset of observed variables and yi,m denote the subset of variables for which realizations are not observable [25, 26]. Finally, ym is the space associated with ym. Without a closed form solution, the maximization of ℓo(θY|Yo=yo) can be cumbersome.

The EM algorithm provides an indirect approach for obtaining the MLE for θY. The process begins by obtaining an estimate of θY using only the completely observed portion of the data, denoted θY(0). Given this initial estimate of the parameters, the E-step finds the expectation of ℓ(θY|y) given the observed data, yo, and the current estimate of θY, θY(0). The M-step maximizes Eℓ(θY|Y)|yo;θY(0) from the E-step, yielding an updated estimate of the parameter, θY(1). The E- and M-steps are repeated with θY(1) used in place of θY(0).

The algorithm iterates through E- and M-steps; the E-step for the (t+1)th iteration is defined as

where θY(t) is the current parameter estimate obtained from the tth iteration. The M-step, requires that we select the updated parameter, θY(t+1) that maximizes the Q-function (Equation (1)) over the parameter space, ΘY, such that

Under the regularity conditions of Dempster et al. [1] and Wu [27], the likelihood, L(⋅), increases with each iteration

requiring only that Equation (1) be continuous in θY. This is a weak condition that is held in many practical situations, and it holds for the exponential family [2].

The E- and M-steps are repeated until the difference between sequential steps in the likelihood changes by an arbitrarily small amount, δ,

Pragmatically, the stopping criterion

is used, probably motivated by Wu [27] who showed

as t→∞, if the limit points are in a compact, connected subset of all stationary points.

The generalized EM algorithm (GEM) relaxes the M-step condition. It requires that Equation (2) holds, permitting the selection of θY(t+1) such that the Q-function increases. This avoids maximizing over all θY∈ΘY, while ensuring that the likelihood does not decrease after each iteration [1]. The limit point of the sequence of estimated parameters is a stationary point of the likelihood, thus providing a local maximum for θY.

Three relevant variations of the EM algorithm are the Monte Carlo EM (MCEM) [28], expectation-conditional maximization (ECM) of Meng and Rubin [22], and the multi-cycle ECM algorithms [22]. The MCEM algorithm replaces a complicated and analytically intractable expectation in the E-step by approximating the integral with a finite sum,

where m denotes the Monte Carlo sample size. In the (t+1)th E-step, the missing data, Ym, is simulated from the current conditional distribution f(Ym|yo;θY(t)) [28]. For each iteration, the ECM has a single E-step followed by a sequence of conditional maximization steps, CM-steps. The constraint function for the CM-steps is defined such that maximization occurs over the entire parameter space [22]. The multi-cycle ECM extends the ECM algorithm. For each iteration, the multi-cycle ECM permits an E-step before a few select CM-steps. The extra computational effort within each iteration results in larger increases in the likelihood, but does introduce the potential for the algorithm to not be a GEM [22].

2.2 Imperfect data and the motivating problem

For subjects i=1,…,n, let Xi denote imperfect data variables as (Xi∗,Xi). If there is only missing data then Xi∗,Xi=Xi,o,Xi,m so that Xi∗=∅. On the other hand, if there is only measurement error then (Xi∗,Xi)=(Xi,o∗,Xi,o) so that Xi,m∗=Xi,m=∅. Finally, if there is both missing data and measurement error, then (Xi∗,Xi)=(Xi,o∗,Xi,m∗,Xi,o,Xi,m). Realizations of measurement error and missing data have similarities; for example, both have at least one variable for which direct observation may not be possible. For measurement error, this is the true value, xi, and for missing data, this is xi,m.

In what follows, we assume that Xi is only measured with error, (Xi∗,Xi)=(Xi,o∗,Xi,o), and that all the measurement error of all subjects’ follows the same distribution and a classical, unbiased, additive, nondifferential measurement error model. The classical, unbiased, additive, nondifferential measurement error model for the ith subject and a random variable, X, is defined as Xi,o∗=Xi,o+εi where εi∼N(0,θM) for i=1,…,n and where θM parameterizes the distribution of εi [29]. We observe xi,o∗ and not xi,o for each of the i units of interest.

The measurement error is unbiased, E(X∗|x)=x, and non-differential (surrogate) defined by the assumption that X∗|x,y=X∗|x or equivalently X∗ and Y are independent given x [10, 29]. This suggests that a surrogate is the measured substitute for the desired variable, but it gives no information beyond that which would have been given by the true variable. When there is more than one mismeasured variable in the model, then we add the further assumption that the associated errors are independent, that is for mismeasured variables Xj,o∗ and Xk,o∗ we assume that εij⊥εik. A measurement error indicator, CiM∈{0,1} with CiM=0 defining the observation of the true value, may be defined as the probability of measurement error in a manner analogous to the following definition for the probability of missingness.

We let Yi denote the imperfect response and assume that the response is subject to missing data, Yi=(Yi,o,Yi,m). Rubin [30] identified three basic types of missing data mechanisms: Missing completely at random (MCAR), missing at random (MAR) and missing not at random (MNAR). For the ith subject, define the missing data indicator as Ci∈{0,1} with Ci=0 defining observation. MCAR occurs when the probability of not observing a response does not depend either the observed or unobserved data, f(Ci|Yi=yi,θC)=f(Ci|θC), where θC parameterizes the missing data mechanism. MAR occurs when the probability of an observation being missing is independent of the unobserved values given the observed values, f(Ci|Yi=yi,θC)=f(Ci|Yi,o=yi,o,θC), where yi,o is the observed value of yi. An MNAR missing data mechanism occurs when the probability of an observation being missing depends on both the observed and unobserved values, f(Ci|Yi=yi,θC).

Assuming unique parameterization of the conditional distributions [9, 26] while restricting our exposition to classical unbiased nondifferential measurement error, MCAR and MAR missing data mechanisms, the joint distribution may be expressed as

where Zi denote complete variables (i.e. no missing data and no measurement error), and θ={θY,θC,θM,θX,θZ}.

The joint distribution (Equation 5) has been decomposed into the product of four conditional distributions and the marginal distribution of the perfect covariates Z. The outcome model, f(Yi|xi,zi;θY), describes the relationship between the variables of interest, (X,Z) and the outcome. The missing data mechanism, f(Ci|xi,zi;θC) is parameterized by θC. The measurement model, f(Xi∗|xi;θM) is the conditional density of X∗ given the unobserved, true value, X=x. The final conditional and marginal distributions comprise the joint distribution of the covariates, sometimes referred to as the exposure model in epidemiological applications.

The outcome model f(Yi|xi,zi;θY) is of primary interest as is an estimate of θY. Typically, the remaining parameters, {θC,θM,θX,θZ}, are treated as nuisance parameters. If the nuisance parameter space is orthogonal to the parameter space of interest, then an unbiased estimator for θY can be obtained provided the nuisance parameter estimator is root-n consistent [23].

We have a primary interest in the outcome model, f(Yi|xi,zi;θY), and a secondary interest in the missing data mechanism, f(Ci|xi,zi;θC). As such, we are interested in the unbiased estimation of {θY,θC}. Under the specification of unique parameterization, the assumption is that {θY,θC} is orthogonal to {θM,θX,θZ}. The challenge is that we are interested in a parameter that is commonly relegated to the nuisance space. An example of possible application (Section 5) is the use of inverse probability weights (IPW) when there is mismeasurement in at least one variable used to model the weights and the response, and the response may be missing [31]. In this scenario, we wish to obtain “good” estimates of the parameters in the weighting models to ensure that the inverse probability weighted estimators exhibit reduced or no bias.

Motivated by this example, we used a simulation study to examine the finite-sample bias reduction for both {θY,θC}, using known implementations of the EM algorithm. Due to lack of convergence of the algorithm for a standard implementation, we implemented the ECM and the multi-cycle ECM algorithm with Monte Carlo integration. By using a single EM algorithm to obtain estimates, we observed that θY and θC were have problems converging. Having a unified approach resulted in a lack of convergence for the single EM algorithm. Convergence rates, across the three implementations of the EM algorithm, ranged from 0% to 60%, suggesting unstable algorithms. To overcome this, we sought to partition the parameter space orthogonally and implement separate EM algorithms for each partition. By having root-n consistent estimators in an particular subspace, we would have unbiased estimators that converge to a stationary point in the complementary subspace.

4.1 Simulation specification

We consider a finite-sample simulation study using the general model given in Equation (5). As noted in Section 2.2, we are interested in the unbiased estimation of {θY,θC}. For the ith subject, we specify the exposure model (Xi,Zi), i=1,⋯,n, to be independently and normally distributed with mean zero, σX,X=σZ,Z=1, and σX,Z=0.2, where {θX,θZ}={μX,μZ,σX,X,σZ,Z,σX,Z}. The response Yi is binary with Yi=1 denoting the occurrence of the event of interest. The missing data indicator for the response, Ci, denotes missing responses as Ci=1.

The response is modeled as a logistic regression,

where θY=(θ0Y,θ1Y,θ2Y). The missing data mechanism is modeled as

where θC=(θ0C,θ1C,θ2C). We set θ0Y=0.40 and θ0C=0.03 and the remaining parameters were randomly generated and restricted to lie in the range (−3,3) for each dimension (Table 1).

Imperfect data are introduced in the following manner. Let X be measured with error using a classical unbiased nondifferential measurement error model, X∗=X+ε. Here, ε∼N(0,θM). Bias reduction as a function of θM is explored to consider the effect of imprecision. We let θM take values {0.1,0.3,0.5,0.7}. With σX,X=1, θM is directly interpreted as the proportion of imprecision in the measurement of X [9].

In applied settings, θM may be estimated through the collection of auxiliary information or by inferring a measurement error model from prior studies. This approach prevents the problem of having a nonidentifiable model ([9], section 4.1). Given that we are not generating auxiliary information from which estimates of θM may be derived, we specify θM for all simulations.

We let K=2 for our simulation study for a two-stage OPEM. The first stage targets θC with {θM,θX,θZ} as nuisance parameters and θY fixed at its initialization value, θ0Y. This stage is an implementation of the MCEM algorithm, given θ0Y. The second stage estimates θY given the stage one estimates.

For each parameterization, we implemented S=100 simulations and each simulation had n=500 observations. We replicated this at each of the four specifications of measurement error noise, θM. For the Monte Carlo integration, we used a burn-in of 1,000 and a Monte Carlo sample of m=2500 as determined by a small Gibbs sampler convergence study and relevant recommendations [33, 39]. We used an MCEM algorithm for each stage of the 2-stage OPEM, and a dissimilarity metric based on a one step lag (e.g. Equation (3)) with measure less than 0.0025 to terminate the algorithm.

Due to the low probability of convergence using the MCEM, ECM and multi-cycle ECM algorithm with Monte Carlo integration (Section 2.2), we compare the OPEM with the naive analysis that assumes complete data. This is a reasonable contrast since the OPEM is not a direct extension of the ECM or multi-cycle ECM, but a new option for implementing an EM algorithm when convergence is problematic. The proposed OPEM algorithm has the advantage of convergence when alternate implementations of the EM algorithm have a low probability of convergence. All simulations were implemented in R [50].

4.2 Simulation results

For the estimated coefficients θˆ1,naiveC and θˆ1,naiveY of the naive analysis, we expect either attenuation, estimates that are smaller in magnitude than the true value, or augmentation, estimates that are larger in magnitude than the true value ([9], section 2.6). Furthermore, we expect an enhanced effect due to the modest correlation between X and Z ([9], Section 2.2) and induced bias for θˆ2,naiveC and θˆ2,naiveY ([10], section 5.3). If the OPEM algorithm performs as theory suggests in Section 3, we should observe bias reduction for the affected estimators without inducing bias in those that were largely unaffected by measurement error. Finally, we note that all of the simulations converged which is a substantial improvement over the convergence rates for standard implementations of the EM algorithm to this problem (Section 2.2).

Although our two-stage OPEM algorithm first estimated {θˆC,θˆM,θˆX} and then estimated θˆY given {θˆC,θˆM,θˆX}, we first present the results for response model since it is of primary interest. We follow the response model results with those of secondary interest, the parameter estimation results for the missing data mechanism.