Computing Pattern-Specific Likelihoods

From observed outcomes to a likelihood

In a pattern-mixture analysis, the outcome distribution is modeled separately by dropout pattern. Computing a likelihood therefore begins with two bookkeeping decisions: identify each subject’s pattern, and identify which outcome components are actually observed for that subject. The likelihood contribution uses the distribution for that pattern and the subject’s observed components; it does not insert unobserved outcomes as if they had been measured.

The goal here is practical as well as mathematical. A correct formula can still produce a wrong fit if the pattern index is misassigned, a variance is allowed to become negative, or an optimizer reports convergence at a poor solution. We will make the contribution explicit, calculate it in small examples, and then use those calculations to motivate implementation checks.

Definition: Let \(R_i=r\) denote subject \(i\)'s dropout pattern, and let \(O_i\) be the set of outcome occasions observed for that subject. If the pattern-specific outcome model has density \(f_r\), the subject’s observed-outcome likelihood contribution is \(f_{r,O_i}(y_{i,O_i}\mid x_i;\theta_r)\). Here \(f_{r,O_i}\) is the density of the observed components under the pattern-\(r\) model, \(x_i\) denotes included covariates, and \(\theta_r\) denotes that model’s parameters.

For independent subjects, the conditional outcome likelihood is the product of these contributions:

L_{\mathrm{out}}(\theta)=\prod_{i=1}^{N} f_{r_i,O_i}(y_{i,O_i}\mid x_i;\theta_{r_i}), \qquad \ell_{\mathrm{out}}(\theta)=\sum_{i=1}^{N}\log f_{r_i,O_i}(y_{i,O_i}\mid x_i;\theta_{r_i}).

The independence assumption is across subjects, conditional on the modeled covariates and pattern-specific parameters. Within a subject, repeated outcomes may be correlated; their joint density must account for that dependence. The likelihood above is conditional on the observed pattern assignments. If pattern probabilities are also modeled, their likelihood is an additional component, not a replacement for the outcome likelihood.

For a multivariate normal pattern-specific model with mean vector \(\mu_r(x_i)\) and covariance matrix \(\Sigma_r\), let \(\mu_{r,O_i}(x_i)\) be the mean subvector and \(\Sigma_{r,O_i,O_i}\) the covariance submatrix for observed occasions. The contribution is the multivariate normal density with exactly those observed components and that submatrix. This is a marginal density for the observed subvector, not a conditional density given the missing components.

Calculating contributions carefully

A scalar normal model gives a useful implementation check. If an observed response \(Y\) in pattern \(r\) has distribution \(N(\mu_r,\sigma_r^2)\), the log-likelihood for observed values \(y_1,\ldots,y_n\) in that pattern is

\ell_r(\mu_r,\sigma_r^2) =-\frac{n}{2}\log(2\pi)-\frac{n}{2}\log(\sigma_r^2) -\frac{1}{2\sigma_r^2}\sum_{j=1}^{n}(y_j-\mu_r)^2, \qquad \sigma_r^2>0.

The variance constraint is part of the model. For unconstrained numerical optimization, a standard parameterization is \(\eta_r=\log(\sigma_r^2)\), so \(\sigma_r^2=\exp(\eta_r)>0\) for every real \(\eta_r\). This avoids asking an optimizer to search over invalid negative variances.

Worked Example: Pattern-Specific Normal Likelihoods

Suppose a hypothetical study has five subjects and one outcome of interest. Three subjects are in pattern A, with observed values 1, 4, and 7; two are in pattern B, with observed values 2 and 4. Fit a separate normal mean and variance in each pattern by maximum likelihood. The values are illustrative, not real study data.

In pattern A, the mean is \(\hat{\mu}_A=(1+4+7)/3=4\). The maximum-likelihood variance divides the sum of squared deviations by \(n=3\):

\hat{\sigma}_A^2 =\frac{(1-4)^2+(4-4)^2+(7-4)^2}{3} =\frac{9+0+9}{3}=6.

The fitted log-likelihood is therefore \(\ell_A=-\frac{3}{2}\log(2\pi)-\frac{3}{2}\log(6)-\frac{18}{12}\). Using \(\log(2\pi)=1.83788\) and \(\log(6)=1.79176\), this is \(\ell_A=-6.94445\), to five decimal places.

In pattern B, \(\hat{\mu}_B=(2+4)/2=3\), and \(\hat{\sigma}_B^2=((2-3)^2+(4-3)^2)/2=1\). Thus \(\ell_B=-\log(2\pi)-\log(1)-2/2=-2.83788\). The conditional outcome log-likelihood across the two patterns is \(\ell_{\mathrm{out}}=\ell_A+\ell_B=-9.78233\).

These are maximum-likelihood variances, calculated with denominator \(n\), not unbiased sample variances calculated with denominator \(n-1\). The latter may be useful for descriptive summaries, but substituting them into this likelihood calculation would not give its maximum-likelihood estimates.

For longitudinal outcomes, the same bookkeeping is essential, but the density is multivariate. For a \(k\)-dimensional normal observation \(y\), its log density is \(-\frac{k}{2}\log(2\pi)-\frac{1}{2}\log|\Sigma|-\frac{1}{2}(y-\mu)^\mathsf{T}\Sigma^{-1}(y-\mu)\). Both the dimension \(k\) and the covariance matrix must match the observed components.

Worked Example: A Multivariate Contribution for Observed Occasions

Consider a hypothetical pattern-specific model for outcomes at three occasions. For a subject in this pattern, only occasions 1 and 2 are observed, with \(y_{O}=(12,9)^\mathsf{T}\). Suppose the model gives mean \(\mu=(10,12)^\mathsf{T}\) at those occasions and the corresponding covariance submatrix is \(\Sigma_O=\begin{pmatrix}4&1\\1&9\end{pmatrix}\). The third occasion is not included in this subject’s contribution.

The determinant is \(|\Sigma_O|=4(9)-1(1)=35\), and the inverse is \(\Sigma_O^{-1}=\frac{1}{35}\begin{pmatrix}9&-1\\-1&4\end{pmatrix}\). With \(y_O-\mu=(2,-3)^\mathsf{T}\), the quadratic form is \((2,-3)\Sigma_O^{-1}(2,-3)^\mathsf{T}=2.4\).

The log-density contribution is \(-\log(2\pi)-\frac{1}{2}\log(35)-\frac{1}{2}(2.4) =-1.83788-1.77767-1.20000=-4.81555\). Its density is approximately \(\exp(-4.81555)=0.00810\). This calculation uses a two-dimensional normal density because two outcomes were observed. Using the full three-occasion density, or conditioning on the unobserved third outcome without a specified value, would not represent this observed-data contribution.

Optimization and constraints

In a real fit, parameters may include pattern-specific regression coefficients, variances, and covariance terms. Some are unconstrained; others must satisfy restrictions such as positive variances and positive-definite covariance matrices. A convenient covariance parameterization is \(\Sigma_r=L_rL_r^\mathsf{T}\), where \(L_r\) is lower triangular and its diagonal entries are constrained positive, often by representing them as exponentials of unconstrained parameters. This guarantees positive definiteness when the diagonal entries are strictly positive.

For a scalar normal model, write \(\eta=\log(\sigma^2)\). For observations \(y_j\), the log-likelihood becomes \(-\frac{n}{2}\log(2\pi)-\frac{n}{2}\eta-\frac{1}{2}\exp(-\eta)\sum_j(y_j-\mu)^2\). Differentiating with respect to \(\mu\) and \(\eta\) gives \(\partial\ell/\partial\mu=\sum_j(y_j-\mu)/\sigma^2\) and \(\partial\ell/\partial\eta=-n/2+\sum_j(y_j-\mu)^2/(2\sigma^2)\). At an interior maximum with nonzero residual variation, setting the first derivative to zero gives \(\hat{\mu}=\bar{y}\); setting the second to zero then gives \(\hat{\sigma}^2=n^{-1}\sum_j(y_j-\bar{y})^2\). This derivation justifies the maximum-likelihood calculation used above.

Worked Example: Checking a Numerical Maximum

Use pattern A from the first example. At \(\mu=4\) and \(\sigma^2=6\), the residuals are \(-3,0,3\), which sum to zero, and their squared sum is 18. The two score components are therefore \(\partial\ell/\partial\mu=0/6=0\) and \(\partial\ell/\partial\eta=-3/2+18/(2\cdot6)=0\). Both first-order checks agree with the fitted values.

Now consider the alternative candidate \(\mu=5,\sigma^2=7\). The squared residual sum is \((1-5)^2+(4-5)^2+(7-5)^2=16+1+4=21\). Its log-likelihood is \(-\frac{3}{2}\log(2\pi)-\frac{3}{2}\log(7)-21/14=-7.17568\). This is lower than the fitted value \(-6.94445\), as expected for this candidate. The comparison is a useful check, although comparing two points alone does not prove that an optimizer has found the global maximum.

A practical fitting routine should be checked in more than one way. Recalculate the log-likelihood from the returned parameters using the same observed components and pattern assignments, but an independently written evaluation if feasible. Use several starting values; compare the final log-likelihoods and fitted parameters. Check that constraints hold, the likelihood is finite, and the gradient is small in the optimizer’s parameterization. A small gradient is evidence of a stationary point, not by itself proof of a global maximum.

Assemble the data by pattern.
For each subject, verify the pattern label, observed occasion set, and covariate values used in the model.

Evaluate stable contributions.
Compute log densities rather than multiplying many small densities. Use linear solves and log determinants rather than explicitly forming a matrix inverse where numerical software permits.

Constrain and optimize.
Use transformations or constrained parameterizations for positive variances and valid covariance matrices; optimize the negative log-likelihood if the software minimizes.

Check the result.
Compare multiple starts, inspect convergence and gradients, and independently recompute the final objective and constraint values.

Pattern frequencies, alternatives, and limitations

If the analysis also models the probabilities of the dropout patterns, let \(\pi_r=P(R_i=r)\), with \(\sum_r\pi_r=1\). Under independent subjects, the joint likelihood factors into the pattern-probability component and the conditional outcome component:

L_{\mathrm{joint}}(\pi,\theta) =\prod_{i=1}^{N}\pi_{r_i}\, \prod_{i=1}^{N} f_{r_i,O_i}(y_{i,O_i}\mid x_i;\theta_{r_i}).

For the five subjects in the first example, if pattern probabilities are estimated from the observed pattern counts, their maximum-likelihood estimates are \(\hat{\pi}_A=3/5=0.6\) and \(\hat{\pi}_B=2/5=0.4\). The pattern-frequency log-likelihood is \(3\log(0.6)+2\log(0.4)=-3.36506\). Under this factorization, the joint log-likelihood is \(-9.78233-3.36506=-13.14739\). The pattern-frequency term must be included only when that component is part of the specified likelihood; reporting it separately prevents confusion about what was optimized.

Direct likelihood maximization differs from weighted estimating equations, which target estimating equations rather than a fully specified probability density. It also differs from Bayesian estimation, which combines a likelihood with a prior and summarizes a posterior distribution. These approaches can share pattern-specific outcome models, but their objective functions and inferential interpretations are not interchangeable. The likelihood calculation itself does not resolve what can be inferred about unobserved outcomes in a dropout pattern; that depends on the model and identifying assumptions.

Common failure modes include using the wrong pattern’s parameter vector, including unobserved responses in a contribution, using a full covariance matrix when only a subvector is observed, and mixing a sample variance with a maximum-likelihood variance. Covariance matrices that are not positive definite can make determinants or log densities invalid. Sparse patterns may also provide little information for estimating a large number of pattern-specific parameters; an optimizer’s convergence message cannot create information absent from the data.

A further limitation is that a high likelihood does not establish that the chosen distribution fits well. Numerical optimization checks ask whether the stated objective was evaluated and optimized reliably. Assessing distributional fit and residual behavior is a separate task. Similarly, uncertainty for pattern-weighted summaries is not obtained merely by reporting the maximized likelihood.

Key takeaway: A pattern-specific likelihood is built from each subject’s observed outcomes under that subject’s pattern model. Correct indexing, valid parameter constraints, and independent numerical checks are essential; pattern-frequency terms belong in the likelihood only when they are explicitly modeled.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions below.

A subject has observed outcomes at occasions 1 and 3, but not occasion 2. Which covariance matrix enters a multivariate normal contribution?
Why is \(\eta=\log(\sigma^2)\) a useful optimization parameterization?
For pattern counts \(n_A=4\) and \(n_B=6\), what are the maximum-likelihood pattern probabilities if the pattern probabilities are modeled as unrestricted multinomial probabilities?
Does a small gradient alone show that a numerical optimizer found the global maximum? Explain.

Solutions. (1) Use the \(2\times2\) covariance submatrix for occasions 1 and 3, and the matching two-element mean subvector; the missing occasion is not inserted into this marginal contribution. (2) Since \(\sigma^2=\exp(\eta)>0\) for every real \(\eta\), the transformation enforces positive variance without a boundary search over negative values. (3) The estimates are \(\hat{\pi}_A=4/10=0.4\) and \(\hat{\pi}_B=6/10=0.6\). (4) No. A small gradient supports the conclusion that the result is near a stationary point, but multiple starts and objective checks are also needed, and a stationary point need not be a global maximum.