2.1 A latent class model for neuropsychological test scores

Suppose we have n IID observations, indexed by $i=1,\ldots ,n$ . Let $\textbf {X}_i = (X_{i,1}, X_{i,2}, \ldots , X_{i,p})^\top \in \mathcal {X}\subset \mathbb {R}^p$ be the covariates and $\textbf {Y}_i = (Y_{i,1}, Y_{i,2},\ldots , Y_{i,d})^\top \in \mathcal {Y} \subset \mathbb {Z}^d$ be the outcome variables of interest, where $\mathcal {Y}$ is a bounded discrete set. In the present study, $\textbf {Y}_i$ represents the neuropsychological test scores of individual i for a total of d total tests, but this framework can be applied to other problems as well. In the article, for notational convenience, we may elect to drop the index i to simplify notation and refer to the generic random variable that is not associated with a specific individual. The jth test score/outcome variable $Y_j$ belongs to the set $\{0,1,2,\ldots ,N_j\}$ . We assume that the covariates are always observed and that only the outcomes may be subject to missingness.

We start with a simpler scenario in which the data are completely observed (no absence of any outcome variables). The more complicated setting with missing outcome variables is addressed in Section 3. When we have multiple discrete outcome variables, the nonparametric approach is not feasible due to the curse of dimensionality. To see this, suppose we have four outcome variables, and each outcome variable has ten possible values. The joint distribution has a dimension of $10^4$ , surpassing the order of magnitude of moderately sized samples. This motivates to use a parametric model.

To model the dependency among discrete variables, we consider a mixture model, which is based on the idea of the LCM. Originally presented in its modern form by Goodman (Reference Goodman1974), the LCM takes the form of the following mixture of multivariate multinomials distribution:

(2.1) $$ \begin{align} p(y; \theta) = \sum_{k=1}^K w_k p_k(y; \theta_k) = \sum_{k=1}^K w_k \prod_{j=1}^d p_{k,j}(y_{j}; \theta_{k,j}), \end{align} $$

where $p_{k,j}$ denotes a multinomial with $\theta _{k,j}$ being its the parameter and the $w_k$ ’s are probability weights for each component such that $\sum _k w_k = 1$ and $w_k>0$ for all k. As $p_{k,j}$ is a multinomial distribution, $\theta _{k,j}$ refers to the vector of probabilities (one for each level of $Y_j$ ) that sum to 1. There is a popular implementation for fitting these models through the poLCA package (Linzer & Lewis, Reference Linzer and Lewis2011).

Let Z denote the component assignment such that $Z = \ell $ if and only if the associated observation is generated from mixture component $\ell $ . Namely, $P(Y=y|Z= k;\theta ) = \prod _{j=1}^d p_{k,j}(y_{j}; \theta _{k,j})$ . By construction, the LCM implicitly assumes that for $j\neq j'$ , the variables $Y_j$ and $Y_{j'}$ are conditionally independent given the latent class Z. This is known as the local independence assumption and is commonly used in the literature to model dependence between discrete random variables. While they are conditionally independent, the variables $Y_j$ and $Y_{j'}$ are allowed to depend marginally for $j\neq j'$ . The local independence assumption makes it convenient to calculate the conditional distributions, which are used in Section 3 when we discuss imputation. Some effort has been made to relax the local independence assumption, such as mixtures of log-linear models (Bock, Reference Bock, Bozdogan, Sclove, Gupta, Haughton, Kitagawa, Ozaki and Tanabe1994) and hierarchical LCMs (Chen et al., Reference Chen, Zhang, Liu, Poon and Wang2012; Zhang, Reference Zhang2004), but these have not been as widely adopted due to the increase in the number of parameters.

Throughout this article, we will refer to the data set containing tuples of the form $(\textbf {X},\textbf {Y})$ as the latent incomplete (LI) data and the data set containing tuples of the form $(\textbf {X},\textbf {Y},Z)$ as the latent complete (LC) data. In the setting of no missing data, the LI data are what we typically have access to, but this is often insufficient for direct model fitting due to the unobserved latent variable Z. We will show that one can estimate the model parameters using an EM algorithm with the LI data. We use this nomenclature to avoid confusion in Section 3 when we encounter the traditional type of missingness with the outcome variables.

There has been previous work to incorporate covariates in the LCM (Huang & Bandeen-Roche, Reference Huang and Bandeen-Roche2004; Vermunt, Reference Vermunt2010). This is related to the mixture of experts models that arose from the machine learning literature (Jacobs et al., Reference Jacobs, Jordan, Nowlan and Hinton1991; Jordan & Jacobs, Reference Jordan and Jacobs1993; Yuksel et al., Reference Yuksel, Wilson and Gader2012). These models generalize the classical mixture model formulation by allowing the model parameters $w_k$ and $\theta _k$ to depend on covariates. The simple mixture of experts model (Gormley & Frühwirth-Schnatter, Reference Gormley and Frühwirth-Schnatter2018) takes the following form:

$$ \begin{align*} p(y | x; \boldsymbol{\theta}) = \sum_{k=1}^K w_k(x;\beta_k) p_k(y; \theta_k), \end{align*} $$

where $\beta _k$ is introduced to allow the weights to depend on covariates. As written, the covariates only adjust the weights placed on each mixture component but do not affect the parameters of the individual component distribution.

To apply these ideas to the neuropsychological data, we combine the mixture of experts and LCM, leading to the following mixture of binomial product experts:

(2.2) $$ \begin{align} p(y|x; \boldsymbol{\beta}, \boldsymbol{\theta}) = \sum_{k=1}^K w_k(x;\beta_k) \underbrace{\prod_{j=1}^d \underbrace{\binom{N_j}{y_j} (\theta_{k,j})^{y_j} (1-\theta_{k,j})^{N_j-y_j}}_{p_{k,j}(y_j;\theta_{k,j})}}_{p_k(y;\theta_k)}, \end{align} $$

where $w_k(x;\beta _k) = \exp (\beta _k^\top \cdot (1,x)) / (\sum _{k'\in [K]} \exp (\beta _{k'}^\top \cdot (1,x)))$ is from the multiple-class logistic regression model. Like in the classical LCM, we utilize the local independence assumption as each mixture decomposes as a product of binomials. However, there is a critical difference between our model and the classical LCM. We use the binomial distribution instead of the multinomial distribution, and this leads to a large reduction in the number of parameters. Since test scores are ordered variables, we expect the reduction to a binomial distribution to be reasonable as scores further away from the mean score can be less likely. Similar to the simple mixture of experts model, each mixture weight depends on covariates.

Here, we will treat $k=1$ as the reference group, so $\beta _1$ is the vector of all $0$ by assumption and does not need to be estimated. For $k \in \{2,3,\ldots ,K\}$ , $\beta _k \in \mathbb {R}^{p+1}$ , and for $k\in [K]$ , $\theta _k \in \mathbb {R}^d$ , we assume that the $w_k$ utilize the logistic function, but it is possible to use other link functions, such as the probit and other monotonic functions (Ouyang & Xu, Reference Ouyang and Xu2022). In our work, the outcome variables are neuropsychological test scores, and we treat the latent groups as subgroups of the population of varying cognitive ability. The five levels of the CDR score offer a natural choice for the number of mixture components.

In our model construction, we view each mixture component as representing a subgroup of the population, where the subgroups have varying cognitive ability that is summarized by $\boldsymbol {\theta }$ . Therefore, we interpret the $\boldsymbol {\theta }$ ’s as attributes of population-level groups and do not allow them to depend on covariates. Since the weights $w_k(x;\beta _k)$ depend on covariates, every individual has their own weights associated with each population-level group. These covariates can be interpreted as a way to construct weights on an individual level for each of these population-level subgroups. A side product of this assumption is that this reduces model complexity. If the $\boldsymbol {\theta }$ ’s depend on covariates, then the number of parameters increases from $O(K(d + p))$ to $O(Kdp)$ , which is a dramatic change even for moderately sized d and p and thus another reason we do not allow this as part of our model. Because the outcome variables are test scores, we can also interpret the quantity $N_j \cdot \theta _{k,j}$ as the mean of the jth test score of the kth latent group, and this is useful as a summary statistic for a given latent class.

Binomial product distributions have previously arisen in the literature. One such application was for modeling test scores for spatial tasks in child development in Thomas et al. (Reference Thomas, Lohaus and Brainerd1993), but they considered a fairly limited setting with no covariates, two components, and two outcome variables. Binomial product and Poisson product distributions have also been used in the statistical ecology literature to model species abundance (Brintz et al., Reference Brintz, Fuentes and Madsen2018; Haines, Reference Haines2016; Kéry et al., Reference Kéry, Royle and Schmid2005). These settings differ from ours because they often treat the total number counts $N_j$ as a random quantity. In item response theory (IRT), the Rasch model posits that the probability of answering each question correctly depends on its difficulty and the individual’s ability, assuming a Bernoulli distribution for each question (Rasch, Reference Rasch1960). However, in our data, only total test scores are available. Thus, assuming uniform difficulty across test questions, we model the scores using Binomial distributions. To our knowledge, our article is the first time binomial products has been used to analyze neuropsychological test score data while incorporating covariates.

Note that for the kth component, $\beta _k \in \mathbb {R}^{p+1}$ and $\theta _k \in \mathbb {R}^d$ , there are $(K-1)(p+1)$ parameters associated with the covariates $\boldsymbol {\beta }$ and $Kd$ parameters associated with $\boldsymbol {\theta }$ and a total of $K(p+d+1)-p-1$ parameters for the whole model. Thus, for fixed covariate and outcome dimensions, the number of parameters grows linearly in K.

2.2 Model fitting

We now describe our procedure for fitting the model. An intuitive approach is to estimate the parameters via the maximum likelihood approach. However, computing the maximum likelihood estimator (MLE) is a nontrivial task because the LI log-likelihood is not concave (Lemma 1 in the Supplementary Material). As such, solving the first-order conditions is no longer sufficient for determining the global maximizer.

Due to Lemma 1 (stated formally in Appendix H.1 of the Supplementary Material), maximizing the LI log-likelihood function directly is not straightforward. Here is an interesting insight from the mixture model literature: if we had access to the latent class label, then the MLE would be easy to construct. This insight motivates us to develop an EM algorithm (Dempster et al., Reference Dempster, Laird and Rubin1977) by augmenting the observed data with the true component label Z to calculate the MLE for $(\boldsymbol {\beta }, \boldsymbol {\theta })$ . For each $i=1,\ldots ,n$ , let $Z_i = \ell $ if the ith observation comes from mixture component $\ell $ . Such data augmentation allows one to bypass the problem of taking the logarithm of a sum. To see this, let $\textbf {X}_{1:n} = (\textbf {X}_1, \ldots , \textbf {X}_n)$ and $\textbf {Y}_{1:n} = (\textbf {Y}_1, \ldots , \textbf {Y}_n)$ be the covariates and outcome variables of all n individuals, respectively. The LC log-likelihood function is written as

(2.3) $$ \begin{align} \ell_{\text{LC},n}&(\beta, \theta; \textbf{X}_{1:n}, \textbf{Y}_{1:n}, \textbf{Z}_{1:n}) = \log \left( \prod_{i=1}^n \prod_{k=1}^K \left(w_k(\textbf{X}_i; \beta_k) p_k(\textbf{Y}_i; \theta_k)\right)^{1(Z_i=k)}\right) \nonumber \\ &= \sum_{i=1}^n\sum_{k=1}^K 1(Z_i=k) \log w_k(\textbf{X}_i;\beta_k) + \sum_{i=1}^n\sum_{k=1}^K 1(Z_i=k)\log p_k(\textbf{Y}_i;\theta_k). \end{align} $$

In the EM algorithm, we start from an initial guess $(\hat {\boldsymbol {\beta }}^{(0)}, \hat {\boldsymbol {\theta }}^{(0)})$ and iterate between an expectation step (E-step) and a maximization step (M-step) to update our guess until convergence. Let t denote the tth iteration of the EM algorithm. In the E-step, we compute the expected value of the complete log-likelihood $\ell _{\text {LC}}(\boldsymbol {\beta }, \boldsymbol {\theta }; \textbf {X},\textbf {Y},\textbf {Z})$ conditional on the observed data $\textbf {X}$ and $\textbf {Y}$ . The expected value of the LC log-likelihood given the observed data forms the Q function in the EM algorithm

$$\begin{align*}Q_{\text{LC}}(\boldsymbol{\beta}, \boldsymbol{\theta} | \textbf{X}, \textbf{Y}; \hat{\boldsymbol{\beta}}^{(t)}, \hat{\boldsymbol{\theta}}^{(t)}) := \mathbb{E}[\ell_{\text{LC}}(\boldsymbol{\beta}, \boldsymbol{\theta}; \textbf{X},\textbf{Y},Z) | \textbf{X}, \textbf{Y}; \hat{\boldsymbol{\beta}}^{(t)}, \hat{\boldsymbol{\theta}}^{(t)}] . \end{align*}$$

In practice, we do not have access to the true expectation, so we work with the sample analog $Q_{\text {LC},n}$ . The sample analog can be expressed as

(2.4) $$ \begin{align} &Q_{\text{LC}, n}(\boldsymbol{\beta}, \boldsymbol{\theta} | \textbf{X}_{1:n}, \textbf{Y}_{1:n}; \hat{\boldsymbol{\beta}}^{(t)}, \hat{\boldsymbol{\theta}}^{(t)}) \nonumber \\ &:= \frac{1}{n}\sum_{i=1}^n\sum_{k=1}^K \hat{\pi}_k^{(t)}(\textbf{X}_i, \textbf{Y}_i) \log w_k(\textbf{X}_i;\beta_k) + \frac{1}{n}\sum_{i=1}^n\sum_{k=1}^K \hat{\pi}_k^{(t)}(\textbf{X}_i, \textbf{Y}_i) \log p_k(\textbf{Y}_i;\theta_k), \end{align} $$

where

(2.5) $$ \begin{align} \hat{\pi}_k^{(t)}(\textbf{X}_i, \textbf{Y}_i) := P(Z=k|\textbf{X}_i,\textbf{Y}_i; \hat{\boldsymbol{\beta}}^{(t)}, \hat{\boldsymbol{\theta}}^{(t)}) = \frac{w_k(\textbf{X}_i;\hat{\beta}_k^{(t)}) p_k(\textbf{Y}_i; \hat{\theta}_k^{(t)})}{\sum_{k'=1}^K w_{k'}(\textbf{X}_i;\hat{\beta}_{k'}^{(t)}) p_{k'}(\textbf{Y}_i; \hat{\theta}_{k'}^{(t)})} \end{align} $$

is the weight of observation $(\textbf {X}_i, \textbf {Y}_i)$ for the kth mixture component.

Since the LC log-likelihood function (2.3) decomposes as the sum of a function of $\boldsymbol {\beta }$ and a function of $\boldsymbol {\theta }$ , the maximization step of the EM algorithm is separable. Note that since $w_k$ and $p_k$ are logistic regression and binomial product models, respectively, one can apply standard off-the-shelf model fitting procedures after reweighting each observation i by $\hat {\pi }_k^{(t)}(\textbf {X}_i,\textbf {Y}_i)$ . In the maximization step, the new estimate for $\boldsymbol {\theta }$ is updated as follows for each $k\in [K]$ and $j \in [d]$ ,

(2.6) $$ \begin{align} \hat{\theta}^{(t+1)}_{k, j} := \dfrac{\sum_{i=1}^n \frac{Y_{i,j}}{N_j} \hat{\pi}_k^{(t)}(\textbf{X}_i, \textbf{Y}_i)}{\sum_{i=1}^n \hat{\pi}_k^{(t)}(\textbf{X}_i, \textbf{Y}_i)}, \end{align} $$

which is simply the standard MLE formed by the sample proportion, reweighted by $\hat {\pi }_k^{(t)}(\textbf {X}_i,\textbf {Y}_i)$ . We update our estimate of $\beta _k$ by fitting a multi-class logistic regression by reweighting each corresponding observation by $\hat {\pi }_k^{(t)}(\textbf {X}_i, \textbf {Y}_i)$ . This is equivalent to regressing the variable $W^{(t)}$ on the covariates X, where $W_i^{(t)} = (\hat {\pi }_1^{(t)}(\textbf {X}_i, \textbf {Y}_i), \ldots , \hat {\pi }_K^{(t)}(\textbf {X}_i, \textbf {Y}_i)) \in S^{K-1}$ for all i, where $S^{K-1}$ is the $(K-1)$ th probability simplex. Note that this is not the standard logistic regression because the outcome variable belongs to a probability simplex. This logistic regression can still be fit using gradient descent in standard R packages.

Algorithm 1 describes the process of fitting the mixture of binomial product experts model in the presence of completely observed covariate and outcome data. For a single initialization point, the EM algorithm is not guaranteed to converge to the global maximizer. In practice, we run Algorithm 1 many times with multiple random initial guesses to explore the parameter space sufficiently and choose the parameter estimate from all of them that maximizes the log-likelihood.

2.3 Identifiability

We now make some brief comments on the identifiability of our model. Identifiability means that each parameter value corresponds to a unique probability distribution. That is, the mapping from the parameter space to the space of probability distributions is one-to-one. If identifiability does not hold, estimation becomes problematic because a unique MLE may not exist. However, this notion of identifiability may be too strong for many practical purposes. For example, even the common problem of label swapping violates this identifiability definition, so researchers often consider the idea of identifiability up to permutation of the parameters. We now consider the notion of generic identifiability, which relaxes the definition of identifiability even further. Generic identifiability means identifiability holds almost everywhere in the parameter space (Allman et al., Reference Allman, Matias and Rhodes2009). More formally, this means that the mapping from the parameter space to the space of probability distributions may fail to be one-to-one only on a set of Lebesgue measure zero. Generic identifiability (holding almost everywhere) can be viewed as an intermediate assumption between two common assumptions: global identifiability (holding everywhere) and local identifiability (holding in a neighborhood of the true parameter). In a practical sense, this means that any such model fit on a given data set is unlikely to be unidentifiable, and we consider this notion to be sufficient for our applied data setting.

We note that the sufficient conditions outlined in Proposition 1 are fairly mild. Intuitively, one would not want the number of parameters to be too large when the data dimension d is small. The bound in Assumption A2 ensures that the number of mixture components K remains appropriately controlled relative to d and $N_j$ , as a large K increases model complexity. When there are at least $d=3$ outcome variables of moderate range, K can still be fairly large. Assumption A3 is to ensure identifiability of the logistic regression parameters. We provide further comments on the assumptions and the proof of this proposition in Appendix B of the Supplementary Material. These assumptions are fairly mild, and we show that they can be met using relatively straightforward generating processes through our simulations in Section 6.

3.1 Missing at random and an imputation strategy

Rubin (Reference Rubin1976) outlined three types of missingness mechanisms: missing completely at random (MCAR), MAR, and MNAR. MCAR assumes that the missingness of the variable is independent of all variables in the data. MAR assumes that the missingness of a variable can only depend on observed variables. MNAR assumes that the missingness of a variable can depend on the value of the variable subject to missingness.

In the MCAR data, model fitting is straightforward because one can run Algorithm 1 on the complete cases, and the estimates of the parameters will be consistent. However, MCAR assumes that the missingness is irrelevant to observed outcomes, which is clearly false in the NACC data as a possible reason to miss some test scores is that the individual was too sick to finish the test. Therefore, we consider the MAR assumption. Formally, the definition of MAR is as follows.

Definition 1 (Missing at random)

The outcome variables $\textbf {Y}$ are MAR if the probability of missingness is dependent only on the variables that are observed for the given pattern. This assumption is written as

(3.1) $$ \begin{align} P(\textbf{R}=r | \textbf{X}, \textbf{Y}) \stackrel{\text{MAR}}{=} P(\textbf{R}=r | \textbf{X}, \textbf{Y}_r) \end{align} $$

for all $r\in \mathcal {R}$ .

Notice that the left-hand side of (3.1) represents the selection probability $P(\textbf {R}=r|\textbf {X},\textbf {Y})$ . This quantity is strictly unidentifiable because it depends on unobserved data, specifically $\textbf {Y}_{\bar {r}}$ . Thus, it cannot be estimated even given infinite observed data. The MAR assumption makes $P(\textbf {R}=r|\textbf {X},\textbf {Y})$ identifiable because it equates it to a probability $P(\textbf {R}=r|\textbf {X},\textbf {Y}_r)$ that can be estimated from the observed data; the variables under the conditioning are all of the variables strictly observed under pattern r. The definition of MAR implies that the probability of a given missing pattern does not depend on variables that are unobserved under that pattern. A key advantage of this assumption is that we avoid the challenge of modeling the joint distribution between $\textbf {Y}$ and $\textbf {R}$ since the missingness mechanism $P(\textbf {R}=r|\textbf {X},\textbf {Y})$ does not need to be modeled directly; this is known as the ignorability property (see the discussion later). Thus, we do not have to deal with making potentially unreasonable modeling assumptions on either the selection model $P(\textbf {R}=r|x,y)$ or the extrapolation distribution $p(y_{\bar {r}} | y_r, x, \textbf {R}=r)$ .

The MAR assumption is untestable because its validity relies strictly on data that is unobserved, and so it cannot be rejected by the observed data. Our primary reason for selecting this assumption is for modeling. This may lead to easier interpretability for scientists and practitioners since we fit a global model $p(y|x)$ across all missing patterns rather than a local model $p(y|x, R = r)$ for every pattern r. In the Alzheimer’s disease literature, having more statistically sound work can be meaningful because of limitations in existing work that we have described previously in Section 1.3.1. Since we believe that it may be plausible because these tests are correlated due to underlying cognitive ability, we can use it as a starting point. This can be viewed as a baseline before proceeding with more complex MNAR assumptions. We recognize that MNAR may be more reasonable since missing test scores may be attributed to sickness. However, this is in itself a very open research question because MNAR is a broad class of assumptions, and performing mixture modeling with MNAR is not straightforward; we leave this for future work.

An additional feature of the MAR assumption is that this assumption offers a simple approach to impute the data, which makes the computation of the MLE a lot easier. Before describing the procedure of updating model parameters, we first introduce a multiple imputation procedure in Algorithm 2 that can be used to fill in the missing data. In the algorithm, for notational convenience, when we write a binary vector in the summation or product, this means we iterate over all indices whose elements are nonzero. For instance, when we write “For j in $(1,0,0,1)$ ,” this is equivalent to “For $j=1,4$ .” Therefore, if $\textbf {R}_i = (1,0,0,1)$ , then “For j in $\textbf {R}_i$ ” corresponds to “For $j=1,4$ ” as well.

As stated, Algorithm 2 describes how to impute the outcome variables, assuming that a good estimate $(\hat {\boldsymbol {\beta }}, \hat {\boldsymbol {\theta }})$ is available. We discuss how to actually obtain such an estimate in Section 3.2. The multiple imputation algorithm exploits the fact that the conditional distribution $p(y_{\bar {r}}|y_r,x)$ for any $r\in \mathcal {R}$ remains a mixture model. There are two steps to performing multiple imputation: 1) for every observation i with missing observations, we compute the weights of each component of the mixture distribution $p(y_{\bar {\textbf {R}}_i}|\textbf {Y}_{i,\textbf {R}_i},\textbf {X}_i)$ and 2) sample M times from the distribution $p(y_{\bar {\textbf {R}}_i}|\textbf {Y}_{i,\textbf {R}_i},\textbf {X}_i)$ for each observation i to form M completed data sets. The derivation of this procedure is provided in Appendix A.2 of the Supplementary Material.

3.2 Model fitting under a missing at random assumption

For each r, we assume that the selection probability $P(\textbf {R}=r|\textbf {X},\textbf {Y}_r; \gamma _r)$ belongs to a parametric family, indexed by $\gamma _r$ . We collect these parameters into $\boldsymbol {\gamma } = (\gamma _r)_{r\in \mathcal {R}}$ . For simplicity, we write the log-likelihood in terms of the probability model without the n samples. Under the MAR assumption, we can write the observed log-likelihood as

$$ \begin{align*} \ell_{\text{obs}}(\boldsymbol{\gamma}, \boldsymbol{\beta}, \boldsymbol{\theta}; x, y_r, r) &= \log P(\textbf{R}=r | x,y_r; \gamma) + \log p(y_r | x; \boldsymbol{\beta}, \boldsymbol{\theta}) \\ &= \ell_{\text{obs}}^{(1)}(\boldsymbol{\gamma}; x, y_r, r) + \ell_{\text{obs}}^{(2)}(\boldsymbol{\beta}, \boldsymbol{\theta}; y_r, x). \end{align*} $$

The missingness mechanism is said to be ignorable because estimation of $(\boldsymbol {\beta }, \boldsymbol {\theta })$ is separated from the estimation of $\boldsymbol {\gamma }$ . Model fitting under an MAR assumption typically occurs using an EM algorithm approach. We can proceed with estimating $(\boldsymbol {\beta }, \boldsymbol {\theta })$ by conditioning on the observed data $(\textbf {X},\textbf {Y}_{\textbf {R}})$ and using another EM algorithm approach. The population $Q_{\text {FC}, r}$ function writes as follows:

(3.2) $$ \begin{align} Q_{\text{FC}, r}(\boldsymbol{\beta}, \boldsymbol{\theta} | \textbf{X}, \textbf{Y}_r; \hat{\boldsymbol{\beta}}^{(t)}, \hat{\boldsymbol{\theta}}^{(t)}) &:= \mathbb{E}[\ell_{\text{LI}}(\boldsymbol{\beta}, \boldsymbol{\theta}; \textbf{X}, \textbf{Y}) | \textbf{X}, \textbf{Y}_r; \hat{\boldsymbol{\beta}}^{(t)}, \hat{\boldsymbol{\theta}}^{(t)}] \nonumber \\ &= \sum_{y_{\bar{r}}} \ell_{\text{LI}}(\boldsymbol{\beta}, \boldsymbol{\theta}; \textbf{X}, y_{\bar{r}}, \textbf{Y}_r) p (y_{\bar{r}} | \textbf{Y}_r, \textbf{X}; \hat{\boldsymbol{\beta}}^{(t)}, \hat{\boldsymbol{\theta}}^{(t)}) \end{align} $$

for every $r \in \mathcal {R}$ . One major observation is that the conditional expectation relies on being able to fit estimate $p(y_{\bar {r}} | y_r, x; \boldsymbol {\beta }, \boldsymbol {\theta })$ . We can leverage the consequences of the local independence assumption to obtain an easily computable form of the conditional distribution. Unfortunately, since the LI log-likelihood function $\ell _{\text {LI}}$ is not linear in the outcome variables, we are unable to reduce (3.2) to a more simple form. However, we can approximate this expectation stochastically by imputing the missing data $\textbf {Y}_{\bar {r}}$ for every missing pattern r using the distribution $p(y_{\bar {r}} | y_r, x; \hat {\boldsymbol {\beta }}^{(t)}, \hat {\boldsymbol {\theta }}^{(t)})$ .

To overcome estimating the conditional expectation of factorial terms, we approximate (3.2) stochastically with a Monte Carlo procedure. For large enough M, we expect that

$$ \begin{align*} Q_{\text{FC}, r}(\boldsymbol{\beta}, \boldsymbol{\theta} | \textbf{X}, \textbf{Y}_r; \hat{\boldsymbol{\beta}}^{(t)}, \hat{\boldsymbol{\theta}}^{(t)}) &\approx Q_{\text{FC},r,n}^{(M)}(\boldsymbol{\beta}, \boldsymbol{\theta} | \textbf{X}_{1:n}, \textbf{Y}_{1:n,\textbf{R}_{1:n}}; \hat{\boldsymbol{\beta}}^{(t)}, \hat{\boldsymbol{\theta}}^{(t)}) \\ &:= \frac{1}{Mn}\sum_{m=1}^M \sum_{i=1}^n \log \left( \sum_{k=1}^K w_k(\textbf{X}_i; \beta_k) p_k(\tilde{\textbf{Y}}_{i}^{(m;t)}; \theta_k)\right) \cdot 1(\textbf{R}_i = r), \end{align*} $$

where $\tilde {\textbf {Y}}_i^{(m;t)} = (\textbf {Y}_{i,\textbf {R}_i}, \tilde {\textbf {Y}}_{i,\bar {\textbf {R}}_i}^{(m;t)})$ denotes the mth imputed data for the ith observation on the iteration step t given the observed variables. The choice of the number of imputations M is important, and in practice, we use $M=20$ to balance both computational time and good estimation performance. We discuss this in more detail with our simulations in Section 6 and Appendix E.2 of the Supplementary Material.

In the maximization step, we compute the MLE of $(\boldsymbol {\beta }, \boldsymbol {\theta })$ using the completed data after the multiple imputation. Since we have completed data, the MLE can be found using the original EM algorithm, outlined in Algorithm 1. We summarize the procedure in the following algorithm. Throughout this article, we will use the notation $\hat {\boldsymbol {\beta }}$ and $\hat {\boldsymbol {\theta }}$ (without superscripts relating to t) to denote the point estimate obtained by Algorithm 3 after convergence is achieved. In practice, we assume that convergence is achieved if the $L_2$ norm of the difference between the current and old parameter estimates falls within a prespecified tolerance level (we choose $\epsilon =10^{-4}$ ).

Algorithm 3 is a nested EM procedure because we have two types of missingness: missingness in the outcome variables and missingness in the latent class labels. We embed an EM algorithm for the LCM fitting (this serves as the M-step in the outer EM algorithm) within an outer EM algorithm for handling the MAR data. Because we do not obtain a closed-form expression for the conditional expectations in the E-step but rather perform a stochastic approximation, Algorithm 3 is a Monte Carlo EM algorithm (Levine & Casella, Reference Levine and Casella2001). We use the notation $\tilde {\textbf {Y}}_i^{(m;t)}$ to denote the mth imputed outcome variables using the model parameterized by $(\hat {\boldsymbol {\beta }}^{(t)}, \hat {\boldsymbol {\theta }}^{(t)})$ . When the parameter is understood, we omit the t index. Note that when there is no missingness in the outcome variables, Algorithm 3 reduces to Algorithm 1 because the multiple imputation step is bypassed. Again, in practice, we run Algorithm 3 with multiple random initial points to ensure we explore the parameter space and converge to the global maximizer. We pick the point estimate that maximizes the observed log-likelihood, defined as

(3.3) $$ \begin{align} \ell_{\text{obs},n} := \ell_{\text{obs},n}(\boldsymbol{\beta}, \boldsymbol{\theta}; \textbf{X}_{1:n}, \textbf{Y}_{1:n, \textbf{R}_{1:n}}) := \sum_{i=1}^n \log \left( \sum_{k=1}^K w_k(\textbf{X}_i) p_{k,\textbf{R}_i}(\textbf{Y}_{i,\textbf{R}_i})\right), \end{align} $$

where $p_{k,r}(y_r) = \prod _{j\in r} p_{k,j}(y_k)$ for any $r\in \mathcal {R}$ . Note that when there is no missing data, this reduces to the LI log-likelihood function $\ell _{\text {LI},n}$ . Since we are working under a parametric model and estimating parameters using maximum likelihood, the estimators converge at $\sqrt {n}$ -rate and are asymptotically efficient. We provide theoretical justification for these methods and discussion of the asymptotic behavior in Appendix A of the Supplementary Material.

7.1 Description of outcome variables and covariates

Our primary goal is to measure the cognitive ability of Alzheimer’s disease patients while measuring its association with baseline covariates. Following similar analysis in Brenowitz et al. (Reference Brenowitz, Hubbard, Keene, Hawes, Longstreth, Woltjer and Kukull2017), we focus on four cognitive domains: long-term memory (episodic memory), attention, language, and executive function. Specifically, we used the Logical Memory Story A immediate and delayed recall to assess memory, Digit Span Forward and Backward tests to evaluate attention, Animal listing test to measure language ability, and Trail Making Tests Parts A and B to measure executive function. Additionally, we used the mini-mental state examination (MMSE) to assess overall cognitive impairment, resulting in a total of eight outcome variables. The Trail A and B tests are assessed based on the time it takes to complete a given task, so a higher score is considered worse. On the other hand, for the other tests, higher scores are indicative of stronger cognitive ability.

We include four baseline covariates in this study: age, education, sex, and race. Age was kept on a yearly scale, and education was dichotomized based on whether an individual had obtained a college degree. Of particular interest was the association between education and cognitive decline. This analysis was motivated by the cognitive reserve hypothesis, which suggests that some mechanisms can provide individuals more resilience to cognitive decline. We would like to explore this hypothesis to see if individuals with higher levels of education are more resistant to cognitive impairment (Meng & D’Arcy, Reference Meng and D’Arcy2012; Thow et al., Reference Thow, Summers, Saunders, Summers, Ritchie and Vickers2018).

For each individual, we have a Clinical Dementia Rating (CDR) score, which is assigned by a medical professional. The CDR score is based on a clinician’s judgment on an individual’s cognitive ability and is the standard procedure to determine whether someone has dementia or not. A CDR score 0 refers to cognitively normal, a score 0.5 refers to an MCI state, and a score of at least 1 refers to dementia. Note that the CDR score is all based on a clinician’s judgment, which is very different from the neuropsychological tests (exam-based). This score takes values in the set $\{0,0.5,1,2,3\}$ . Here, 0 indicates normal cognitive ability, 0.5 is MCI, 1 is mild dementia, 2 is moderate dementia, and 3 is severe dementia. Thus, we expect that a healthy group would contain individuals with lower CDR scores than an unhealthier group. We will use the CDR score to interpret some of our results, but we emphasize that we do not use this variable in the model.

Since there are eight outcome variables, there are up to $2^8 = 256$ missing data patterns. In this data set, we observe 93 missing data patterns with seven patterns comprising over 95% of the data. In particular, there are $n_{cc}=25,041$ fully observed cases and $n=41,181$ total observations, so that approximately only 60% of the individuals have completely observed outcomes.

In Figure 2, we plot the CDR distribution for the complete cases, the individuals missing at least one outcome variable, and the entire data set, respectively. If the data were MCAR, we would expect the CDR score distribution to be mostly constant across all missing patterns of the data. However, we observe that for the complete cases, the dementia group (CDR score of 1, 2, or 3) makes up less than 20%, but the dementia group is close to 40% of the data that has at least one missing outcome variable. This suggests that the dementia group is severely underrepresented when we only perform analysis on the complete cases. Thus, from the plots, we can visually postulate the MCAR assumption to be unreasonable.

Figure 2 The CDR score distributions for the complete cases, the individuals missing at least one outcome variable, and the entire data set are provided in the left, middle, and right panels, respectively.

We select the model using prior knowledge. Since the scientific community previously determined that five levels for the CDR score were appropriate, we also choose the same number. A similar number of groups will help with interpretability, and five is small enough that the number of parameters is manageable.

7.2 Point estimates and confidence intervals

The data analysis was conducted on a 2022 Macbook Air with 8 cores, and it took about a day to fit the model and complete the bootstrap with 1,000 bootstrap samples. In Table 4, we report the point estimates of the test score means and the corresponding standard errors in parentheses. There are a few takeaways. As reported, Classes 1 and 5 can be interpreted as the most and least healthy individuals, respectively. Note that in 6 out of 8 tests (omitting TRAIL A and TRAIL B), a higher score corresponds to better performance. On the other hand, TRAIL A and TRAIL B are scored using time-to-completion in seconds, so a higher score corresponds to a worse performance. For every test aside from TRAIL B, we observe a clear monotonic behavior in the test score means from Class 1 to Class 5, which suggests that the model is reasonable and fairly interpretable. For TRAIL B, Classes 4 and 5 have very similar mean test scores that are close to the maximum possible test score of 300. Since a test score of 300 indicates that the individual timed out and did not complete the test, this suggests that individuals from Classes 4 and 5 have comparable performance on the TRAIL B test.

Table 4 These tables contain point estimates of the test score means for each latent class

Note: In the title of each column, we report the maximum possible score for a given test in parentheses. In each cell, the standard errors are reported in parentheses.

Additionally, the results also suggest the different tests may have value at distinguishing individuals from different latent cognitive classes. For example, TRAIL A and TRAIL B are routinely used to measure executive function. Classes 4 and 5 represent the most unhealthy individuals in the population, but individuals in both perform similarly on TRAIL B. On the other hand, the mean test scores in Classes 4 and 5 for TRAIL A are very different, which suggests that TRAIL A can be a good discriminator between very unhealthy and moderately unhealthy individuals. As TRAIL B is a more complex task than TRAIL A, this seems to agree with our intuition.

In Table 5, we report the point estimates of the coefficients in the logistic regressions. One initial observation is that the coefficients of Age and Education are all positive and negative at the $\alpha =0.05$ significance level, respectively, across all classes. This implies that increasing age is associated with lower cognitive ability and higher education may provide some protection against dementia. The latter is known as the cognitive reserve hypothesis in the context of education (Meng & D’Arcy, Reference Meng and D’Arcy2012; Thow et al., Reference Thow, Summers, Saunders, Summers, Ritchie and Vickers2018). We can compare the magnitudes of the coefficient of Education to the coefficient of Age. For example, in Class 2, we observe $0.62/0.06 \approx 10$ , and in Class 5, we observe $1.14/0.05\approx 23$ . This suggests that having a college degree or higher education may be equivalent to being 10–20 years younger in terms of cognitive ability.

Table 5 This table contains point estimates for the coefficients for each of the covariates and each latent class

Note: In each cell, the standard errors are reported in parentheses.

7.3 Latent classes and clustering

We report two clustering results using models fit on the complete cases (Figure 3) and the entire data (Figure 4). For each set of clustering results, we visualize the clusters by plotting the CDR score composition of each of the five clusters via barplots. Furthermore, we report the mean of the CDR scores for each group, which gives a rough estimate on the group-specific clinical cognitive ability. From the complete case results (CC) in Figure 3, we generally see that the first clusters are primarily composed of healthy individuals as individuals with CDR scores of 0 take up larger proportions. As we progress to the middle latent classes, healthier individuals take up fewer proportions and individuals with MCI (CDR score of 0.5) take up larger proportions. Then, in Group 5, the majority of the individuals have some form of dementia.

Figure 3 Clustering on complete data only.

Note: These barplots summarize the composition of each of the five latent groups. We order the groups from the most healthy to the least healthy. This is reflected in the mean CDR score of each group.

Figure 4 Clustering on the entire data with MAR assumption.

Note: These barplots summarize the composition of each of the five latent groups. We order the groups from the most healthy to the least healthy. This is reflected in the mean CDR score of each group.

For the clustering that used the entire data with the MAR assumption (Figure 4), we observe a similar trend with the healthier individuals, the MCI individuals, and the dementia individuals being the dominant proportion in the earlier, middle, and later latent classes. However, the results from clustering the entire data set are more stark. In contrast with the complete clustering, Groups 4 and 5 contain a significantly higher number of individuals. Moreover, our model is able to detect a group of individuals with high dementia, as evidenced by Group 5’s mean CDR score of 1.55. Unlike the complete case clustering, Group 5 has almost 75% of individuals with some form of dementia. Many of these individuals would be omitted from the analysis if the missing data were not accounted for. Thus, accounting for the missingness leads to a clustering result with a stronger correlation with clinical assessments (CDR score).