Worked Gaussian Missing-Data Analysis

From Missing Values to a Gaussian Analysis

A small data set with several missingness patterns is useful for seeing what a Gaussian likelihood actually uses. A complete record contributes a joint density; a record with only one component observed contributes that component’s marginal density. In EM, the unobserved components enter through their conditional moments, not through arbitrary single-value substitutions.

We work with a bivariate normal model \(X_i\sim N_2(\boldsymbol\mu,\boldsymbol\Sigma)\), where \(\boldsymbol\Sigma\) is positive definite. The observed-data likelihood below assumes the missingness mechanism is ignorable for likelihood inference: the missingness is MAR and its parameters are distinct from the parameters of the Gaussian outcome model. The analysis is about that outcome-model likelihood; it does not establish that this assumption is appropriate in a particular study.

Definition: For record \(i\), let \(O_i\) denote the set of observed components and \(M_i\) the set of missing components. The observed-data likelihood contribution is the marginal density of \(X_{i,O_i}\) under the model. In an EM iteration, the E-step computes conditional expectations of complete-data sufficient statistics given the observed values and current parameter values; the M-step maximizes the resulting expected complete-data log-likelihood.

Conditional Moments for an Incomplete Gaussian Record

Partition a Gaussian vector into an observed block \(X_O\) and a missing block \(X_M\), with corresponding partitions of its mean and covariance. When \(\boldsymbol\Sigma_{OO}\) is nonsingular, the conditional distribution of the missing block is

X_M\mid X_O=x_O \sim N\left( \boldsymbol\mu_M+\boldsymbol\Sigma_{MO}\boldsymbol\Sigma_{OO}^{-1}(x_O-\boldsymbol\mu_O), \boldsymbol\Sigma_{MM}-\boldsymbol\Sigma_{MO}\boldsymbol\Sigma_{OO}^{-1}\boldsymbol\Sigma_{OM} \right).

For a single missing scalar \(Y\) and observed scalar \(X=x\), write its conditional mean and variance as \(m_{Y\mid x}\) and \(v_{Y\mid x}\). Then the conditional moments needed for a bivariate Gaussian E-step include

E[Y\mid X=x]=m_{Y\mid x},\qquad E[Y^2\mid X=x]=v_{Y\mid x}+m_{Y\mid x}^2,\qquad E[XY\mid X=x]=x\,m_{Y\mid x}.

The last identity holds because \(X=x\) is fixed after conditioning. The second is especially important: replacing \(Y\) by its conditional mean alone discards its conditional variance and therefore understates the expected second moment.

Result and derivation: The stated conditional distribution follows by defining \(R=X_M-\boldsymbol\mu_M-\boldsymbol\Sigma_{MO}\boldsymbol\Sigma_{OO}^{-1}(X_O-\boldsymbol\mu_O)\). Its covariance with \(X_O\) is \(\boldsymbol\Sigma_{MO}-\boldsymbol\Sigma_{MO}\boldsymbol\Sigma_{OO}^{-1}\boldsymbol\Sigma_{OO}=0\). Because \(R\) and \(X_O\) are jointly Gaussian, zero cross-covariance implies independence. Direct covariance expansion gives \(\operatorname{Var}(R)=\boldsymbol\Sigma_{MM}-\boldsymbol\Sigma_{MO}\boldsymbol\Sigma_{OO}^{-1}\boldsymbol\Sigma_{OM}\). Thus conditioning on \(X_O=x_O\) fixes the affine term and leaves the independent residual \(R\), proving the formula. \(\square\)

Worked Example: A Conditional Second Moment

Suppose \(\boldsymbol\mu=(10,5)^\mathsf{T}\) and \(\boldsymbol\Sigma=\begin{pmatrix}4&1.2\\1.2&2.25\end{pmatrix}\), and the first component is observed as \(X=13\), while \(Y\) is missing. Compute the conditional mean, variance, and second moment of \(Y\).

The conditional mean is \(5+(1.2/4)(13-10)=5+0.9=5.9\). The conditional variance is \(2.25-(1.2^2/4)=2.25-0.36=1.89\). Therefore \(E[Y^2\mid X=13]=1.89+5.9^2=1.89+34.81=36.70\). The expected cross-product is \(E[XY\mid X=13]=13(5.9)=76.7\). These are the conditional moments required for the corresponding sufficient-statistic entries.

A Small Data Set and Two EM Updates

Consider five records, with the two components denoted \(X\) and \(Y\):

Data: Three complete records are \((1,0)\), \((0,2)\), and \((2,1)\). A fourth record has \(X=1\) and missing \(Y\); a fifth has missing \(X\) and \(Y=1\).

There are eight observed scalar values: six in the complete records and one in each incomplete record. Start EM at \(\boldsymbol\mu^{(0)}=(0,0)^\mathsf{T}\) and \(\boldsymbol\Sigma^{(0)}=\mathbf I_2\). With this starting model, either missing component has conditional mean zero and variance one, whatever the observed component. The conditional expected sufficient statistics for the two incomplete records are consequently

\begin{array}{c|ccccc} & E[X]&E[Y]&E[X^2]&E[Y^2]&E[XY]\\ X=1,\;Y\text{ missing}&1&0&1&1&0\\ X\text{ missing},\;Y=1&0&1&1&1&0 \end{array}

For the complete records, the sums of \(X,Y,X^2,Y^2,XY\) are respectively \(3,3,5,5,2\). Adding the two conditional rows gives expected totals \(4,4,7,7,2\). The Gaussian M-step uses the expected complete-data means and covariance, with divisor \(n=5\):

\boldsymbol\mu^{(1)}= \begin{pmatrix}4/5\\4/5\end{pmatrix},\qquad \boldsymbol\Sigma^{(1)} = \begin{pmatrix}7/5&2/5\\2/5&7/5\end{pmatrix} - \begin{pmatrix}4/5\\4/5\end{pmatrix} \begin{pmatrix}4/5&4/5\end{pmatrix} = \begin{pmatrix}0.76&-0.24\\-0.24&0.76\end{pmatrix}.

This covariance is positive definite: its eigenvalues are \(0.76-0.24=0.52\) and \(0.76+0.24=1.00\). For a second E-step, under these updated parameters, the missing \(Y\) given \(X=1\) has mean \(0.8+(-0.24/0.76)(0.2)=14/19\) and variance \(0.76-0.24^2/0.76=13/19\). By symmetry, the same conditional mean and variance apply to missing \(X\) given \(Y=1\). For either incomplete record, the missing component’s expected square is \(13/19+(14/19)^2=443/361\), and the expected cross-product is \(14/19\).

The complete-record totals stay unchanged. The updated expected totals are therefore \(\sum E[X]=\sum E[Y]=3+1+14/19=90/19\), \(\sum E[X^2]=\sum E[Y^2]=5+1+443/361=2609/361\), and \(\sum E[XY]=2+28/19=66/19\). Applying the same M-step gives

\boldsymbol\mu^{(2)}= \begin{pmatrix}18/19\\18/19\end{pmatrix},\qquad \boldsymbol\Sigma^{(2)} = \begin{pmatrix}989/1805&-366/1805\\-366/1805&989/1805\end{pmatrix}.

This second covariance is also positive definite: its diagonal is positive and its determinant is \((989^2-366^2)/1805^2=844165/1805^2>0\).

Worked Example: Computing the Observed-Data Likelihood

For an observed component set \(O_i\), the marginal distribution is Gaussian with mean \(\boldsymbol\mu_{O_i}\) and covariance \(\boldsymbol\Sigma_{O_iO_i}\). Thus, for nonempty observed sets, the observed log-likelihood is

\ell_{\mathrm{obs}}(\boldsymbol\mu,\boldsymbol\Sigma) = -\frac{1}{2}\sum_{i=1}^n \left[ |O_i|\log(2\pi)+\log|\boldsymbol\Sigma_{O_iO_i}| +(x_{i,O_i}-\boldsymbol\mu_{O_i})^\mathsf{T} \boldsymbol\Sigma_{O_iO_i}^{-1} (x_{i,O_i}-\boldsymbol\mu_{O_i}) \right].

This expression follows by multiplying the Gaussian marginal densities for the observed components of each record and taking logarithms. Integrating the joint Gaussian density over missing components gives that marginal density; no density factor is contributed for an unobserved value itself.

At the starting parameters, the complete records have squared Mahalanobis distances \(1,4,5\), while the two observed singletons each have squared distance \(1\). All observed covariance determinants equal one. The quadratic terms sum to \(12\), and the eight observed scalars contribute \(-8\log(2\pi)/2=-4\log(2\pi)\). Hence \(\ell_{\mathrm{obs}}(\boldsymbol\mu^{(0)},\boldsymbol\Sigma^{(0)})=-4\log(2\pi)-6\approx-13.35\).

At the first update, the determinant of the complete-record covariance is \(0.76^2-0.24^2=0.52\). For the three complete records, the residual vectors from the updated mean are \((0.2,-0.8)\), \((-0.8,1.2)\), and \((1.2,0.2)\). Using \(\boldsymbol\Sigma^{-1}=0.52^{-1}\begin{pmatrix}0.76&0.24\\0.24&0.76\end{pmatrix}\), their quadratic forms are \(0.44/0.52=11/13\), \(1.12/0.52=28/13\), and \(1.24/0.52=31/13\). Their sum is \(70/13\). Each singleton has marginal variance \(0.76\), residual \(0.2\), and quadratic form \(0.04/0.76=1/19\). The log-likelihood is therefore

\ell_{\mathrm{obs}}(\boldsymbol\mu^{(1)},\boldsymbol\Sigma^{(1)}) =-4\log(2\pi)-\frac32\log(0.52)-\log(0.76) -\frac12\left(\frac{70}{13}+\frac{2}{19}\right) \approx -8.84.

The normalizing constant is based on eight observed scalars, not on five records or on the three complete records alone. In particular, each singleton contributes one scalar Gaussian normalizing constant.

Complete-Case Estimates and What They Do Not Tell You

Worked Example: Comparing with the Complete-Case Fit

Using only the three complete records, the complete-case mean is \(\bar X=\bar Y=1\). Their deviations from that mean are \((0,-1)\), \((-1,1)\), and \((1,0)\). With the Gaussian maximum-likelihood divisor \(3\), the complete-case covariance estimate is

\widehat{\boldsymbol\Sigma}_{\mathrm{CC}} = \begin{pmatrix}2/3&-1/3\\-1/3&2/3\end{pmatrix}, \qquad |\widehat{\boldsymbol\Sigma}_{\mathrm{CC}}|=1/3.

This estimate differs from the first EM estimate \(\boldsymbol\Sigma^{(1)}\), which uses all five records through their observed values and conditional moments. To compare these parameter values on the same observed data, evaluate the observed likelihood at the complete-case estimates. The inverse complete-case covariance is \(\begin{pmatrix}2&1\\1&2\end{pmatrix}\). The three complete-record quadratic forms are each \(2\), for a total of \(6\). Both singleton observations equal their fitted marginal means, so their quadratic terms are zero; their marginal variances are \(2/3\). Thus

\ell_{\mathrm{obs}}(\widehat{\boldsymbol\mu}_{\mathrm{CC}}, \widehat{\boldsymbol\Sigma}_{\mathrm{CC}}) =-4\log(2\pi)-\frac32\log(1/3)-\log(2/3)-3 \approx -8.30.

For this data and these parameter values, the complete-case fit has a higher observed likelihood than the first EM iterate, approximately \(-8.30\) versus \(-8.84\). That does not contradict EM’s likelihood-increase property: that property compares consecutive iterates from the same EM run, not an iterate with every other possible parameter value. Nor does this one-iteration comparison show that the complete-case estimate is the observed-likelihood maximizer. Further EM iterations are needed before assessing the converged likelihood-based fit.

Why the Pattern-Specific Calculation Matters

The likelihood calculation illustrates why missingness patterns should be handled explicitly. A complete bivariate record contributes a determinant for the full \(2\times2\) covariance; a singleton contributes only the variance of its observed component. Counting the wrong number of observed scalar components changes the normalizing constant, while using the full covariance determinant for a singleton changes the likelihood itself.

The E-step has a parallel lesson. Conditional means supply first moments, but the conditional covariance must also enter the expected second moments. In this example, the term \(13/19\) in the second update is part of \(E[Y^2\mid X]\); dropping it would produce a different covariance update. The conditional cross-moment must also be computed under the conditional distribution rather than set to zero merely because a value is missing.

Common pitfall: Treating a missing component as equal to its conditional mean and then proceeding as if that value were observed. The correction is to use the full conditional expectation of each complete-data sufficient statistic: for a missing scalar, include its conditional variance in its expected square, and use the conditional mean in the expected cross-product.

Check Your Understanding

Use the data and notation in this tutorial to check each step of the analysis.

How many observed scalar components contribute Gaussian normalizing constants in the five-record data set, and what is their total log-normalizing contribution?
For a missing \(Y\) with conditional mean \(m\) and variance \(v\), why is \(E[Y^2\mid X]=m^2+v\) rather than \(m^2\)?
At the starting parameters, verify the expected totals of \(X^2\), \(Y^2\), and \(XY\) used in the first M-step.
Why does a singleton observation use a marginal variance rather than the determinant of the full bivariate covariance?
Does the complete-case fit’s observed likelihood exceeding the first EM iterate contradict EM monotonicity? Explain the comparison to which monotonicity applies.