2 Inhomogeneous HMM for multivariate longitudinal data with state-dependent missingness

In the remainder, we will consider HMMs with K states for longitudinal data consisting of sets of time-series (e.g., time-series for different patients) which may differ in length. Let $Y_{1:N,1:T_i} = (Y_{1,1:T_1}, Y_{2,1:T_2},\ldots ,Y_{N,1:T_n})$ denote such a set of N time-series $Y_{i, 1:T_i} = (Y_{i,1}, \ldots , Y_{i,T_i})$ , each of length $T_i$ . Whilst we will focus on univariate responses in the remainder, the results apply directly to D-variate responses, $Y_{i,t} = (Y_{i,t,1},\ldots , Y_{i,t,D})$ , as long as conditional independence (12) holds. We will allow state-transitions to be inhomogeneous (i.e., time-variant) by including covariates $\mathbf {x}_{i,1:T_i}$ on the initial state probabilities and state-transition probabilities. Here, we use multinomial logistic regressions:

$$\begin{align*}p(S_{i,1} = j|{\boldsymbol\theta}_{\text{pr}}, \mathbf{x}_{i,1}) = \frac{\exp ({\boldsymbol\beta}_j^{(\text{pr})} \mathbf{x}_{i,1})}{\sum_{k=1}^K \exp ({\boldsymbol\beta}_k^{(\text{pr})} \mathbf{x}_{i,1})},\end{align*}$$

and

$$\begin{align*}p(S_{i,t+1} = k|S_{i,t} = j, {\boldsymbol\theta}_{\text{tr}}, \mathbf{x}_{i,t}) = \frac{\exp ({\boldsymbol\beta}_{j,k}^{(\text{tr})} \mathbf{x}_{i,t})}{\sum_{l=1}^K \exp ({\boldsymbol\beta}_{j,l}^{(\text{tr})} \mathbf{x}_{i,t})}.\end{align*}$$

Note that for identification, ${\boldsymbol \beta }_k^{(\text {pr})}$ and ${\boldsymbol \beta }_{j,k}^{(\text {tr})}$ should be fixed to 0 for one state $k \in (1, \ldots , K)$ , and that usually, $\mathbf {x}_{i,1:T_i}$ will include a constant term for the intercept.

We allow responses $Y_{i,t}$ to depend on hidden states $S_{i,t}$ and covariates $\mathbf {x}_{i,t}$ . For continuous-valued variates $Y_{i,t}$ , we can for example use linear regressions:

$$\begin{align*}p(Y_{i,t}|{\boldsymbol\theta}_{\text{obs}}, \mathbf{x}_{i,t}, S_t = j) = \mathbf{Normal}({\boldsymbol\beta}_{j}^{(\text{obs})} \mathbf{x}_{i,t}, \sigma_{j}).\end{align*}$$

Finally, missingness $M_{i,t}$ is allowed to depend on the hidden states and covariates. Here, we use logistic regression

$$\begin{align*}p(M_{i,t} = 1|{\boldsymbol\phi}, \mathbf{x}_{i,t}, S_t = j) = \frac{\exp ({\boldsymbol\beta}_{j}^{(\text{mis})} \mathbf{x}_{i,t})}{1 + \exp ({\boldsymbol\beta}_{j}^{(\text{mis})} \mathbf{x}_{i,t})}.\end{align*}$$

2.1 Estimation via expectation–maximization

Then the following, let ${\boldsymbol \theta } = ({\boldsymbol \theta }_{\text {pr}}, {\boldsymbol \theta }_{\text {tr}}, {\boldsymbol \theta }_{\text {obs}}, {\boldsymbol \phi })$ denote all the model parameters, with ${\boldsymbol \theta }_{\text {pr}} = ({\boldsymbol \beta }_1^{(\text {pr})},\ldots ,{\boldsymbol \beta }_K^{(\text {pr})})$ denoting the parameters for the initial state probabilities, ${\boldsymbol \theta }_{\text {tr}} = ({\boldsymbol \beta }_{1,1}^{(\text {tr})},\ldots ,{\boldsymbol \beta }_{K,K}^{(\text {tr})})$ the parameters for the state-transition probabilities, ${\boldsymbol \theta }_{\text {obs}} = ({\boldsymbol \beta }_{1}^{(\text {obs})},\ldots ,{\boldsymbol \beta }_{K}^{(\text {obs})},\sigma _{1},\ldots ,\sigma _{K})$ the parameters for the state-conditional observation densities, and ${\boldsymbol \phi } = ({\boldsymbol \beta }_1^{(\text {mis})},\ldots ,{\boldsymbol \beta }_K^{(\text {mis})})$ the parameters for the state-conditional missingness probabilities.

These parameters can be estimated through the expectation–maximization (EM) algorithm, which in the context of HMM is also known as the Baum–Welch algorithm. The EM algorithm consists of iteratively maximizing the expected joint log-likelihood

$$\begin{align*}Q({\boldsymbol\theta}, {\boldsymbol\theta}') = \mathbb{E}[\log p(y_{1:N,1:T}, m_{1:N, 1:T}, s_{1:N,1:T}|{\boldsymbol\theta}) ]\end{align*}$$

where the expectation is taken with respect to $p(S_{1:N,1:T}|y_{1:N,1:T}, m_{1:N,1:T}, {\boldsymbol \theta }')$ . Note that the expectation is based on initial parameter values ${\boldsymbol \theta }'$ , whilst the joint log likelihood is defined over parameter values ${\boldsymbol {\theta }}$ .

The expected joint log-likelihood can be written as

(14) $$ \begin{align} Q({\boldsymbol\theta}, {\boldsymbol\theta}') &= \sum_{i=1}^N \sum_{j=1}^K \gamma_{i,1}(j) \log p(S_{i,1} = j|{\boldsymbol\theta}_{\text{pr}}, \mathbf{x}_{i,1}) \nonumber\\ &\quad+ \sum_{i=1}^N \sum_{t=2}^{T_i} \sum_{j=1}^K \sum_{k=1}^K \xi_{i,t-1}(j,k) \log p(S_{i,t} = k|S_{i,t-1} = j), {\boldsymbol\theta}_{\text{tr}}, \mathbf{x}_{i,t-1}) \nonumber\\ &\quad+ \sum_{i=1}^N \sum_{t=1}^{T_i} \sum_{j=1}^K \gamma_{i,1}(j) \log p(M_{i,t}=1|S_{i,t} = j, \mathbf{x}_{i,t},{\boldsymbol\phi})\nonumber \\ &\quad+ \sum_{i=1}^N \sum_{t=1}^{T_i} \sum_{j=1}^K \gamma_{i,1}(j) \log p^*(y_{i,t}|S_{i,t} = j, \mathbf{x}_{i,t},{\boldsymbol\theta}_{\text{obs}}) \end{align} $$

where

$$\begin{align*}\gamma_{i,1}(j) \stackrel{\text{def}}{=} p(S_{i,t} = j| m_{i,1:{T_i}}, y_{i,1:{T_i}}, \mathbf{x}_{i,1:{T_i}}, {\boldsymbol\phi},{\boldsymbol\theta})\end{align*}$$

is the posterior probability of state $S_{i,t}$ and

$$\begin{align*}\xi_{i,t}(j,k) \stackrel{\text{def}}{=} p(S_{i,t+1} = k, S_{i,t} = j|m_{i,1:{T_i}}, y_{i,1:{T_i}}, \mathbf{x}_{i,1:{T_i}}, {\boldsymbol\phi},{\boldsymbol\theta}_{\text{tr}}),\end{align*}$$

is the joint posterior probability of states $S_{i,t}$ and $S_{i,t+1}$ . These probabilities can be efficiently computed via the forward-backward algorithm (Rabiner, Reference Rabiner1989). We define the forward-variable

$$\begin{align*}\alpha_{i,t}(j) \stackrel{\text{def}}{=} p(m_{i,1:t}, y_{i,1:t},S_{i,t} = j, {\boldsymbol\theta}'),\end{align*}$$

which can be computed iteratively as

(15) $$ \begin{align} \alpha_{i,1}(j) = p(S_{i,1} = j|{\boldsymbol\theta}_{\text{pr}}',\mathbf{x}_{i,1}) p(m_{i,1} | {\boldsymbol\phi}', \mathbf{x}_{i,1},S_{i,1} = j) p^*(y_{i,1} | \mathbf{x}_{i,1}, {\boldsymbol\theta}_{\text{obs}}', S_{i,1} = j) \end{align} $$

and for $t> 1$ , as

(16) $$ \begin{align} \alpha_{i,t}(j) = \sum_{l=1}^{K} \alpha_{i,t-1}(l) p(S_{i,t} = j|S_{i,t-1}=l, \mathbf{x}_{i,t-1}, {\boldsymbol\theta}_{\text{tr}}') p(m_{i,t} | {\boldsymbol\phi}', \mathbf{x}_{i,t}) p^*(y_{i,t} | \mathbf{x}_{i,t}, {\boldsymbol\theta}_{\text{obs}}'). \end{align} $$

We also define the backward-variable

$$\begin{align*}\beta_{i,t}(j) \stackrel{\text{def}}{=} p (m_{i,(t+1):T_i},y_{i,(t+1):{T_i}}|S_{i,t} = j, \mathbf{x}_{i,(t+1):{T_i}}, {\boldsymbol\theta}_{\text{obs}}'),\end{align*}$$

which is initialized at $t=T_i$ as

(17) $$ \begin{align} \beta_{T_i}(j) = 1, \quad j = 1, \ldots, K \end{align} $$

and then for each time $t=T_i - 1, \ldots , 1$ as

(18) $$ \begin{align} \beta_{i,t}(j) &= \sum_{l = 1}^K p(S_{i,t+1} = l| S_{i,t} = j, \mathbf{x}_{i,t}, {\boldsymbol\theta}_{\text{tr}}') \nonumber\\ &\quad\times p(m_{i,t+1}|S_{i,t+1} = l, {\boldsymbol\phi}', \mathbf{x}_{i,t+1}) p^*(y_{i,t+1}|S_{i,t+1} = l, {\boldsymbol\theta}_{\text{obs}}', \mathbf{x}_{i,t+1}). \end{align} $$

Using the forward and backward variables, we can compute the posterior state probabilities as

(19) $$ \begin{align} \gamma_{i,t}(j) = \frac{\alpha_{i,t}(j) \beta_{i,t}(j)}{\sum_{l=1}^K \alpha_{i,t}(l) \beta_{i,t}(l)} \end{align} $$

and

(20) $$ \begin{align} \xi_{i,t}(j,k) = \frac{\alpha_{i,t}(j) p(S_{i,t+1} = k|S_{i,t} = j, {\boldsymbol\theta}_{\text{tr}}', \mathbf{x}_{i,t})\beta_{i,t+1}(k)}{\sum_{j=1}^K \sum_{k=1}^K \alpha_{i,t}(j) p(S_{i,t+1} = k|S_{i,t} = j, {\boldsymbol\theta}_{\text{tr}}', \mathbf{x}_{i,t})\beta_{i,t+1}(k)}. \end{align} $$

Note that the expected joint log-likelihood (14) is the sum of four weighted log-likelihoods, one for each set of parameters ${\boldsymbol \theta }_{\text {pr}}$ , ${\boldsymbol \theta }_{\text {tr}}$ , ${\boldsymbol \theta }_{\text {obs}}$ , and ${\boldsymbol \phi }$ . Maximizing the expected joint log-likelihood therefore consists of separately maximizing four weighted likelihoods. When the initial states, state transitions, responses and missingness indicators are modeled with generalized linear models and multinomial logistic regression models, as we have done here, we can then rely on the standard maximum likelihood estimation procedures for these models (see McCullagh & Nelder, Reference McCullagh and Nelder1989), using the $\gamma _{i,1}(j)$ and $\xi _{i,t}(j,k)$ values as case-weights (see also Visser & Speekenbrink, Reference Visser and Speekenbrink2022).

The full EM algorithm can be specified as

1. 1. Start with initial parameters ${\boldsymbol \theta }'$ .

2. 2. Do until convergence:

   1. a. For $i=1, \ldots , N$ , $t=1, \ldots , T_i$ , $j, k = 1, \ldots , K$ , compute $\gamma _{i,t}(j)$ (19) and $\xi _{i,t}(j,k)$ (20) via the forward-backward recursions (15, 16, 17, 18).

   2. b. Obtain new estimates

      $$\begin{align*}\hat{{\boldsymbol\theta}}_{\text{pr}} = \arg \max_{{\boldsymbol\theta}_{\text{pr}}} \sum_{i=1}^N \sum_{j=1}^K \gamma_{i,1}(j) \log(p(S_{i,1} = j|{\boldsymbol\theta}_{\text{pr}}, \mathbf{x}_{i,1})),\end{align*}$$
      $$\begin{align*}\hat{{\boldsymbol\theta}}_{\text{tr}} = \arg \max_{{\boldsymbol\theta}_{\text{tr}}} \sum_{i=1}^N \sum_{t=2}^{T_i} \sum_{j=1}^K \sum_{k=1}^K \xi_{i,t-1}(j,k) \log p(S_{i,t} = k|S_{i,t-1} = j), {\boldsymbol\theta}_{\text{tr}}, \mathbf{x}_{i,t-1}),\end{align*}$$
      $$\begin{align*}\hat{{\boldsymbol\theta}}_{\text{obs}} = \arg \max_{{\boldsymbol\theta}_{\text{obs}}} \sum_{i=1}^N \sum_{t=1}^{T_i} \sum_{j=1}^K \gamma_{i,1}(j) \log p^*(y_{i,t}|S_{i,t} = j, \mathbf{x}_{i,t},{\boldsymbol\theta}_{\text{obs}}),\end{align*}$$

      and

      $$\begin{align*}\hat{{\boldsymbol\phi}} = \arg \max_{{\boldsymbol\theta}_{\text{tr}}} \sum_{i=1}^N \sum_{t=1}^{T_i} \sum_{j=1}^K \gamma_{i,1}(j) \log p(m_{i,t}|S_{i,t} = j, \mathbf{x}_{i,t},{\boldsymbol\phi}).\end{align*}$$

   3. c. If $\frac {L(\hat {{\boldsymbol \theta }}|Y) - L({\boldsymbol \theta }'|Y)}{L({\boldsymbol \theta }'|Y)} < \epsilon $ (e.g., $\epsilon = 1 \times 10^{-8}$ ), assume convergence.

   4. d. Set ${\boldsymbol \theta }' = (\hat {{\boldsymbol \theta }}_{\text {pr}}, \hat {{\boldsymbol \theta }}_{\text {tr}}, \hat {{\boldsymbol \theta }}_{\text {obs}}, \hat {{\boldsymbol \phi }}).$

The EM algorithm is guaranteed to converge to a local maximum of the likelihood. Assessing whether the algorithm converged to the global maximum is not possible in general. To increase the chances to obtain the global maximum likelihood parameters, the algorithm can be run many times, each time using different starting values ${\boldsymbol \theta }'$ . Starting values for ${\boldsymbol \theta }_{\text {pr}}$ and ${\boldsymbol \theta }_{\text {tr}}$ can be derived by assuming uniform distributions for $p(S_1)$ and $p(S_{t},S_{t-1})$ , or sampling these from suitable Dirichlet distributions. Starting values for ${\boldsymbol \theta }_{\text {obs}}$ are less straightforward to choose in general, and arguably more important. One method is to randomly sample state probabilities $\gamma _{i,t}(j)$ from a suitable Dirichlet distribution, and then set ${\boldsymbol \theta }_{\text {obs}}'$ to the maximum likelihood estimates as in step b. This usually provides valid starting values and is the default option in depmixS4 (Visser & Speekenbrink, Reference Visser and Speekenbrink2010).

2.2 Model selection, checking, and standard errors

An important consideration when using HMMs is the number of latent states K. This is generally determined by estimating models with different values for K and then choosing the best one via model selection criteria such as the Akaike information criterion (AIC; Akaike, Reference Akaike1998) and the Bayesian information criterion (BIC; Schwarz, Reference Schwarz1978). For present purposes, another consideration is choosing between a MAR and MNAR model. As our MNAR model contains and additional missingness variable M, the AIC and BIC measures cannot be used directly, as the models target different likelihoods (one for the joint distribution of Y and M, and one for the distribution of just Y). A suitable alternative is to fix the probability of missingness in the MNAR model to be identical over the states, effectively creating a MAR model. As this MAR model is nested within the MNAR model, a general likelihood ratio test

$$\begin{align*}-2 \log \frac{L({\boldsymbol\theta}_{\text{MNAR}}|Y_{1:N,1:T},M_{1:N,1:T})}{L({\boldsymbol\theta}_{\text{MAR}}|Y_{1:N,1:T},M_{1:N,1:T})} \sim \chi^2\left(|{\boldsymbol\theta}_{\text{MNAR}}| - |{\boldsymbol\theta}_{\text{MAR}}|\right)\end{align*}$$

may be used to determine whether the MNAR model outperforms the MAR model ( $|{\boldsymbol \theta }|$ denotes the number of free parameters in ${\boldsymbol {\theta }}$ ).

Another consideration is whether the distributional assumptions for the state-conditional distributions $p(Y_{i,t}|S_t)$ are reasonable. Zucchini et al. (Reference Zucchini, MacDonald and Langrock2017) propose computing “pseudo-residuals” from the cumulative probabilities

$$\begin{align*}p(Y_t \leq y_t|Y_{1:(t-1)},Y_{(t+1):T}), \end{align*}$$

which are converted to the corresponding quantiles of a standard Normal distribution. If the model fits the data, then these quantiles will follow a standard Normal distribution. They show that the cumulative probability can be computed as a weighted sum

$$\begin{align*}\sum_{j=1}^K w_{i,t,j} p(Y_{i,t} \leq y_{i,t}|S_{i,t} = j)\end{align*}$$

with

$$\begin{align*}w_{i,t,j} \propto \begin{cases} p(S_{i,1} = j) \beta_{i,1}(j) & t = 1 \\ \sum_{k=1}^K \alpha_{i,t-1}(k) p(S_{t+1} = j|S_{t}=k) \beta_{i,t}(j) & t> 1 \end{cases}\end{align*}$$

where the weights are normalized such that $\sum _j w_{i,t,j} = 1$ .

For inference on model parameters, standard errors and confidence intervals are of importance. There are several methods to compute (approximate) standard errors for maximum likelihood parameter estimates of HMM (Visser et al., Reference Visser, Raijmakers and Molenaar2000). The standard approach for obtaining standard errors from the Hessian matrix of second derivatives of the log-likelihood function is computationally tricky, as discussed in Visser et al. (Reference Visser, Raijmakers and Molenaar2000), although Lystig and Hughes (Reference Lystig and Hughes2002) provide an elegant solution to overcome this computational challenge. Other methods include likelihood profiles and bootstrapping. Here we use a finite differences approach to estimate the Hessian matrix, which in turn is used to compute confidence intervals for the estimated parameters. Note that the method for finite differences proposed in Visser et al. (Reference Visser, Raijmakers and Molenaar2000) was updated in Visser and Speekenbrink (Reference Visser and Speekenbrink2022) and implemented in depmixS4 (Visser & Speekenbrink, Reference Visser and Speekenbrink2010). This updated finite difference method provides standard error estimates that are as accurate as those provided by bootstrapping methods (which are much more computationally expensive).

3 Simulation study

To assess the potential benefits of including a state-dependent missingness model in a HMM, we conducted a simulation study, focusing on a three-state HMM with a univariate Normal distributed response variableFootnote ³. We simulated four scenario’s. In Simulation 1 and 2, the states are reasonably well-separated with means $\mu _1 = -1$ , $\mu _2 = 0$ , $\mu _3 = 1$ and standard deviations (SD) $\sigma _1 = \sigma _2 = \sigma _3 = 1$ (see Figure 3). Note that there is still considerable overlap in the state-conditional response distributions, as would be expected in many real applications of HMMs. In Simulation 1, missingness was state-dependent (i.e., MNAR), with $p(M_{i,t} = 1|S_{i,t} = 1) = .05$ , $p(M_{i,t} = 1|S_{i,t} = 2) = .25$ , and $p(M_{i,t} = 1|S_{i,t} = 3) = .5$ . In Simulation 2, missingness was independent of the state (MAR), with $p(M_{i,t} = 1|S_{i,t}=i) = p(M_{i,t}) = .25$ . In Simulation 3 and 4 (Figure 3), the states were less well-separated, with means as for Simulation 1 and 2, but SD $\sigma _i = 3$ (see Figure 3). Here, the overlap of the state-conditional response distributions is much higher than in Simulation 1 and 2, and identification of the hidden states will be more difficult. In Simulation 3, missingness was state-dependent (MNAR) in the same manner as Simulation 1, while in Simulation 4, missingness was state-independent (MAR) as for Simulation 2. In all simulations, the initial state probabilities were $\pi _1 = p(S_{i,1} = 1) = .8$ , $\pi _2 = \pi _3 = .1$ , and the state-transition matrix was

$$ \begin{align*} \mathbf{A} = \left[ \begin{matrix} .75 & .125 & .125 \\ .125 & .75 & .125 \\ .125 & .125 & .75 \end{matrix} \right]. \end{align*} $$

In each simulation, we simulated a total of 1000 data sets, each consisting of $N = 100$ replications of a time-series of length $T = 50$ . We denote observations in such replicated time series as $Y_{i,t}$ , with $i=1,\ldots ,N$ and $t = 1, \ldots , T$ . Data was generated according to a 3-state hidden Markov model. For MAR cases, the non-missing observations are distributed as

(21) $$ \begin{align} p(Y_{i,t}|S_{i,t} = j) = \mathbf{Normal}(Y_{i,t}|\mu_j,\sigma_j). \end{align} $$

In the MNAR cases, the missingness variable M and the response variable Y were conditionally independent given the hidden state:

(22) $$ \begin{align} p(Y_{i,t},M_{i,t}|S_{i,t} = j) = \mathbf{Bernouilli}(M_{i,t}|\phi_j) \times \mathbf{Normal}(Y_{i,t}|\mu_j,\sigma_j). \end{align} $$

Data sets were simulated by first generating the hidden state sequences $S_{i,1:T}$ according to the initial state and transition probabilities. Then, the observations $Y_{i,1:T}$ were sampled according to the state-conditional distributions $p(Y_{i,t}|S_{i,t})$ . Finally, random observations were set to missing values according to the missingness distributions $p(M_{i,t}|S_{i,t})$ .

Figure 3

State-conditional response distributions in the simulation studies. In Simulation 1 and 2, states are reasonably well-separated, although there is still considerable overlap of the distributions. In Simulation 3 and 4, states are less well-separated.

We fitted two 3-state HMM to each data-set. In the MAR models, observed responses were assumed to be distributed according to (21), and in the MNAR models, the observed responses and missingness indicators were assumed to be distributed according to (22). Parameters were estimated by maximum likelihood, using the EM algorithm, as implemented in depmixS4 (Visser & Speekenbrink, Reference Visser and Speekenbrink2010). To speed up convergence, starting values were set to the true parameter values. Although such initialization is obviously not possible in real applications, we are interested in the quality of parameter estimates at the global maximum likelihood solution, and setting starting values to the true parameters makes it more likely to arrive at the global maximum. In real applications, one would need to use a sufficient number of randomly generated starting values to find the global maximum.

The results of Simulation 1 (Table 2) show that, when states are relatively well separated, both models provide parameter estimates which are, on average, reasonably close to the true values. Both models have the tendency to estimate the means as more dispersed, and the SD as slightly smaller, then they really are. While wrongly assuming MAR may not lead to overly biased estimates, we see that the MAE for the MNAR model is always smaller than that of the MAR model, reducing the estimation error to as much as 58%. Over all parameters, the relative MAE of the models is 0.77 on average, which shows a clear advantage of the MNAR model. As such, accounting for state-dependent missingness increases the accuracy of the parameter estimates. We next consider recovery of the hidden states, by comparing the true hidden state sequences to the maximum a posteriori (MAP) state sequences determined by the Viterbi algorithm (see Rabiner, Reference Rabiner1989; Visser & Speekenbrink, Reference Visser and Speekenbrink2022). The MAR model recovers 53.13% of the states, while the MNAR model recovers 62.86% of the states. The accuracy in recovering the hidden states is thus higher in the model which correctly accounts for state-dependent missingness. Whilst the performance of neither model may seem overly impressive, we should note that recovering the hidden states is a non-trivial task when the state-conditional response distributions have considerable overlap (see Figure 3) and states do not persist for long periods of time (here, the true self-transitions probabilities are $a_{ii} = .75$ , meaning that states have an average run-length of 4 consecutive time-points). When ignoring time-dependencies and treating the observed data as coming from a bivariate mixture distribution over Y and M, the maximum accuracy in classification would be 50.09% for this data. The theoretical maximum classification accuracy for the hidden Markov model is more difficult to establish, but simulations show that the MNAR model with the true parameters recovers 66.51% of the true states. For the MAR model, the approximate maximum classification accuracy is 58.06%.

Table 2

Results of Simulation 1 (MNAR, low variance)

Note: Values shown are the true value of each parameter, and the mean (mean), standard deviation (SD), and mean absolute error (MAE) of the parameter estimates, for both the MAR and MNAR model. The value of “rel. MAE” is the ratio of the mean absolute error of the MNAR over the MAR model.

The results of Simulation 2 (Table 3) show that when data is in fact MAR, both models provide roughly equally accurate parameter estimates. Whilst the MNAR model does not provide better parameter estimates, including a model component for state-dependent missingness does not seem to bias parameter estimates compared to the MAR model. As can be seen, the state-wise missingness probabilities are, on average, close to the true values of .25. Over all parameters, the relative MAE of the models is 1.003 on average, which shows the models perform equally well. In terms of recovering the hidden states, the MAR model recovers 55.6% of the states, while the MNAR model recovers 55.63% of the states. The somewhat reduced recovery rate of the MNAR model compared to Simulation 1 is likely due to the fact that here, missingness provides no information about the identity of the hidden state. For comparison, the maximum classification accuracy is 42.91% for a mixture model, and approximately 60.45% for the HMM.

Table 3

Results of Simulation 2 (MAR, low variance)

Note: Values shown are the true value of each parameter, and the mean (mean), standard deviation (SD), and mean absolute error (MAE) of the parameter estimates, for both the MAR and MNAR model. The value of “rel. MAE” is the ratio of the mean absolute error of the MNAR over the MAR model.

In Simulation 3 (Table 4) and 4 (Table 5) the states are less well-separated, making accurate parameter estimation more difficult. Here, the tendency to estimate the means as more dispersed and the SD as smaller than they are becomes more pronounced. For both models the estimation error in Simulation 3 (Table 4) is larger than for Simulation 1, but comparing the MAE for both models again shows the substantial benefits of including a state-dependent missingness model. Over all parameters, the relative MAE of the models is 0.658 on average, which shows the MNAR model clearly outperforms the MAR model. In terms of recovering the hidden states, the MAR model recovers 34.97% of the states, whilst the MNAR model recovers 45.27% of the states. As in Simulation 1, the MNAR model performs better in state identification. For both models, performance is lower than in Simulation 1, reflecting the increased difficulty due to increased overlap of the state-conditional response distributions (Figure 3). Indeed, the performance of the MAR model is close to chance (random assignment of states would give an expected accuracy of 33.33%). The maximum classification accuracy is 44.03% for a mixture model, and approximately 54.04% for the MNAR and 41.42% for the MAR HMM.

Table 4

Results of Simulation 3 (MNAR, high variance)

Note: Values shown are the true value of each parameter, and the mean (mean), standard deviation (SD), and mean absolute error (MAE) of the parameter estimates, for both the MAR and MNAR model. The value of “rel. MAE” is the ratio of the mean absolute error of the MNAR over the MAR model.

Table 5

Results of Simulation 4 (MAR, high variance)

Note: Values shown are the true value of each parameter, and the mean (mean), standard deviation (SD), and mean absolute error (MAE) of the parameter estimates, for both the MAR and MNAR model. The value of “rel. MAE” is the ratio of the mean absolute error of the MNAR over the MAR model.

When missingness is ignorable (Simulation 4), like in Simulation 2, inclusion of a state-dependent missingness component in the HMM does not increase bias in parameter estimates. Over all parameters, the relative MAE of the models is 0.987 on average, which shows the models perform roughly equally well. The MAR model recovers 35.51% of the states, whilst the MNAR model recovers 35.5% of the states. For comparison, the maximum accuracy is 36.64% for a mixture model, and 42.51% for the HMM.

Taken together, these simulation results show that if missingness is state-dependent, there is a substantial benefit to including a (relatively simple) model for missingness in the HMM. When missingness is in fact ignorable, including a missingness model is superfluous, but does not bias the results. Hence, there appears to be little risk associated with including a missingness submodel in the HMM.

Four additional simulations were conducted to assess to what extent these results depend on the persistence of states and the homogeneity of state transition probabilities over the states. In these simulations, we used the relatively well-separated states of Simulations 1 and 2. In Simulation 5 and 6, we changed the initial state probabilities to a uniform distribution $\pi _1 = p(S_{i,1} = 1) = \pi _2 = \pi _3 = 1/3$ and the state-transition matrix to

$$ \begin{align*} \mathbf{A} = \left[ \begin{matrix} .5 & .25 & .25 \\ .25 & .5 & .25 \\ .25 & .25 & .5 \end{matrix} \right]. \end{align*} $$

Thus, in these simulations, initial state identification may be more difficult, and the states are (even) less persistent than in Simulations 1 and 2. Full results are provided in the Appendix. When data is MNAR (Table A1), we again find a clear advantage of the MNAR model, with an average relative MAE of 0.886. The MAR model recovers 44.38% of the states, while the MNAR model recovers 51.62% of the states. When data is MAR (Table A2), both models perform roughly equally well, with an average relative MAE of 0.971. The MAR model recovers 46.33% of the states, while the MNAR model recovers 46.34% of the states. In Simulations 7 and 8, we used the same initial state probabilities as in Simulation 1 and 2, but changed the state-transition matrix to

$$ \begin{align*} \mathbf{A} = \left[ \begin{matrix} .8 & .15 & .05 \\ .0375 & .85 & .1125 \\ .025 & .075 & .9 \end{matrix} \right]. \end{align*} $$

As such, the states persist longer than in Simulation 1 and 2, and persistence is furthermore dependent on the state. When data is MNAR (Table A3), we again find a clear advantage of the MNAR model, with an average relative MAE of 0.727. The MAR model recovers 64.52% of the states, while the MNAR model recovers 73.56% of the states. When data is MAR (Table A4), both models perform roughly equally well, with an average relative MAE of 1.003. The MAR model recovers 68.09% of the states, while the MNAR model recovers 68.03% of the states.

In a further simulation, we consider a more traditional case of MNAR data, where missingness depends on the underlying value of the response variable. More specifically, we model the probability of missingness as a function of the true value of the response $Y_{i,t}$ via a logistic regression:

$$\begin{align*}p(M_{i,t} = 1|Y_{i,t}) = \frac{1}{1+\exp(-1 \times (-2 + 2 \times Y_{i,t}))}, \end{align*}$$

keeping the other parameters the same as in Simulation 1. Whilst missingness does not directly depend on the hidden state, because the true response values do depend on the states, the probability of missingness differs between the states, with approximately 6.8%, 22.5%, and 50% expected missing values in states 1, 2, and 3 respectively. As such, although the relation between the underlying true value of the response and missingness is not part of the MNAR model, we would expect the model to indicate state-dependent missingness. The results of this simulation (Table 6) show a clear bias in estimating the state-dependent means and SD: because the higher the value of the response, the higher the probability that value is missing, both state dependent means and SD are underestimated, particularly for state 3 where the probability of missingness is highest. Whilst bias in parameter estimates is evident in both the MAR and MNAR model, the latter performs better on average: over all parameters, the relative MAE of the models is 0.835 on average, which shows a clear advantage of the MNAR model. State recovery seems relatively unaffected by the bias in parameter estimates. The MAR model recovers 49.6% of the states, and the MNAR model recovers 59.15% of the states. These results are close to those of Simulation 1. Thus, for this more traditional form of MNAR data, the (misspecified) MNAR model again outperforms the MAR model, and state recovery seems mostly unaffected by the unavoidable bias in parameter estimates.

Table 6

Results of Simulation 9 (MNAR, related to true value)

Note: Values shown are the true value of each parameter, and the mean (mean), standard deviation (SD), and mean absolute error (MAE) of the parameter estimates, for both the MAR and MNAR model. The value of “rel. MAE” is the ratio of the mean absolute error of the MNAR over the MAR model.

In a final simulation, we assessed the performance of the models when missingness is time-dependent, rather than state-dependent. Attrition is common in longitudinal studies, meaning that the probability of missingness often increases with time. In this simulation, the probability of missingness as a function of time t was modeled through a logistic regression model:

(23) $$ \begin{align} p(M_{i,t} = 1) = \frac{1}{1+\exp(-(0.125 \times t - 5))}. \end{align} $$

Here, the probability of missing data is very small (0.008) at time 1, but increases substantially to (0.777) at time 50. The other parameters were the same as in Simulation 1 and 2. In a model that specifies missingness as state-dependent, but not time-dependent, this could potentially result in biased parameter estimates. For instance, the increased probability of missingness over time may be accounted for by estimating states to have a different probability of missingness, and estimating prior and transition probabilities to allow states with a higher probability of missingness to occur more frequently later in time. In addition to the two HMMs estimated before, we also estimated a HMM with a state- and time-dependent model for missingness:

(24) $$ \begin{align} p(M_{i,t} = 1|S_{i,t} = j) = \frac{1}{1+\exp(-(\beta_{0,j} + \beta_{\text{time},j} \times t))}. \end{align} $$

This model should be able to capture the true pattern of missingness, whilst the MNAR model which only includes state-dependent missingness would not.

The results (Table 7) show that, compared to the MAR model, the MNAR model which misspecifies missingness as state-dependent is inferior, resulting in more biased parameter estimates. Over all parameters, the relative MAE of these two models is 1.353 on average, indicating the MAR model outperforms the MNAR (state) model. To account for the increase in missing values over time, the MNAR (state) model estimates the probability of missingness as highest for state 2, which is estimated to have a mean of close to 0, but an increased SD to incorporate observations from the other two states. To make state 2 more prevalent over time, transition probabilities to state 2 are relatively low from state 1 and 3 (parameters $a_{12}$ and $a_{32}$ respectively), whilst self-transitions ( $a_{22}$ ) are close to 1 (meaning that once in state 2, the hidden state sequence is very likely to remain in that state). The prevalence of state 2 is thus increasing over time, and as this state has a higher probability of missingness, so is the prevalence of missing values. The MNAR (time) model, which allows missingness to depend on both the hidden states and time, performs only slightly worse than the MAR model, with an average relative MAE over all parameters of this model compared to the MAR of 1.042. However, the MNAR (time) model is able to capture the pattern of attrition (increased missing data over time), whilst the MAR model is not. As such, the MNAR (time) model may be deemed preferable to the MAR model, insofar as one is interested in more than modeling the responses Y. In terms of recovering the hidden states, the MAR model recovers 55.67% of the states, and the MNAR (time) model recovers 55.42% of the states. The misspecified MNAR (state) model recovers 50% of the states. The maximum classification accuracy for this data is 42.95% for a mixture model, and approximately 59.91% for the HMM.

Table 7

Results of Simulation 10 (time-dependent missingness, low variance)

Note: Values shown are the true value of each parameter, and the mean (mean), standard deviation (SD), and mean absolute error (MAE) of the parameter estimates, for the MAR, MNAR (state), and MNAR (time) model. The value of “rel. MAE 1” is the ratio of the mean absolute error of the MNAR (state) over the MAR model, and the value of “rel. MAE 2” is the ratio of the mean absolute error of the MNAR (time) over the MAR model. Note that the SDs of $\beta _0,j$ and $\beta _{\text {time},j}$ are relatively high. Whilst the estimates are generally accurate, there are rare outlying estimates which inflate these SDs.

This final simulation shows that when modeling patterns of missing data in HMMs, care should be taken in how this is done. An increase in missing data over time could be due to an underlying higher prevalence of states which result in more missing data, and/or a state-independent increase in missingness over time. In applications where the true reason and pattern of missingness is unknown, it is then advisable to start by allowing for both state- and time-dependent missing data, selecting simpler options when this is warranted by the data.

4 Application to the National Institute of Mental Health Schizophrenia Collaborative Study

In applying HMMs to the National Institute of Mental Health Schizophrenia Collaborative Study, we assume the severity of schizophrenia is characterized by abrupt – rather than continuous—changes. We fitted HMMs in which we either assumed ratings are MAR, or assume ratings are MNAR and allow missingness to be both state- and time-dependent. For each type of model (MAR or MNAR), we fit versions with 2, 3, 4, or 5 states. Both types of model assume imps79, the IMPS Item 79 ratings, follow a Normal distribution, with a state-dependent mean and SD. No additional covariates were included on these means, as the states are intended to capture all the important determinants of illness severity. To model effects of drug, we allow transitions between states, as well as the initial state, to depend on a dummy-coded covariate drug (1 for medication, 0 for placebo). Whilst the initial measurement at week 0 was made before administering the drug, we allow the initial state at week 0 to depend on drug in order to account for any potential pre-existing differences between the conditions. In the MNAR models, the missingness variable is modeled with a logistic regression, using week (between 0 and 5) and the dummy-coded main variable (1 for main measurement occasion, 0 for the other occasions) as predictors, as these were found to be important predictors in the (state-independent) logistic regression analysis reported earlier (Table 1). All models were estimated by maximum likelihood using the EM algorithm implemented in depmixS4 (Visser & Speekenbrink, Reference Visser and Speekenbrink2010). Table 8 contains the goodness-of-fit statistics for all the fitted models.

Table 8

Modeling results for the MAR and MNAR HMM with 2-5 latent states

Note: Note that the likelihood and hence the AIC and BIC values cannot be compared between the MAR and MNAR models, as the latter are based on the additional missingness variable.

For both the MAR and MNAR models, the BIC indicates a three-state model fits best, whilst the AIC indicates a five-state model (or higher) fits best. Favouring simplicity, we follow the BIC scores here, and focus on the three-state models. Considering the absolute fit to the data, the pseudo-residuals of the three-state MAR and MNAR models (Figure 4) are similar. Whilst there are to-be-expected deviations due to the mostly discrete nature of the ratings, the distribution of the pseudo-residuals is close to standard Normal, indicating a satisfactory fit to the data. As such, there is no reason to doubt the assumption that the IMPS ratings follow a state-conditional Normal distribution.

Figure 4

Histograms and QQ plots of the pseudo-residuals for the MAR and MNAR model.

We first consider the parameter estimates of the MAR model. The estimated means and SD for the severity of symptoms are

$$ \begin{align*} {\boldsymbol\mu} = [2.315,4.339, 5.7] \quad \quad {\boldsymbol\sigma} = [0.821,0.619, 0.567]. \end{align*} $$

Hence, the states are ordered, with state 1 being the least severe, and state 3 the most severe. The prior probabilities of the states, for treatment with placebo and medication respectively, are

$$ \begin{align*} {\boldsymbol\pi}_{\text{placebo}} = [0,0.333,0.667] \quad \quad {\boldsymbol\pi}_{\text{medication}} = [0.005,0.307,0.689], \end{align*} $$

and the transition probability matrices (with initial states in rows and subsequent states in columns) are

$$ \begin{align*} \mathbf{A}_{\text{placebo}} = \left[ \begin{matrix} 0.963 & 0.005 & 0.032 \\ 0.118 & 0.878 & 0.004 \\ 0.027 & 0.046 & 0.927 \end{matrix} \right] \quad \quad \mathbf{A}_{\text{medication}} = \left[ \begin{matrix} 1 & 0 & 0 \\ 0.231 & 0.764 & 0.005 \\ 0.073 & 0.307 & 0.62 \end{matrix} \right]. \end{align*} $$

As expected, the initial state probabilities show little difference between the treatments (as the initial measurement was conducted before treatment commenced), but the transition probabilities indicate that for those who received medication, transitions to less severe states are generally more likely, indicating effectiveness of the drugs. This is particularly marked for the most severe state, where the probability of remaining in that state is 0.927 with placebo, but 0.62 with medication. Also note the difference between the transition probabilities for the least severe state: when administered medication, the probability of remaining in the least severe state equals approximately 1, whereas that is not the case for the placebo group.

We next consider the three-state MNAR model. The means and SD for the severity of symptoms are

$$ \begin{align*} {\boldsymbol\mu} = [2.325,4.424, 5.757] \quad \quad {\boldsymbol\sigma} = [0.833,0.669, 0.547] \end{align*} $$

showing the same ordering of states in terms of severity. The prior probabilities for placebo and medication conditions are

$$ \begin{align*} {\boldsymbol\pi}_{\text{placebo}} = [0,0.394,0.606] \quad \quad {\boldsymbol\pi}_{\text{medication}} = [0.004,0.349,0.647], \end{align*} $$

and the transition probability matrices are

$$ \begin{align*} \mathbf{A}_{\text{placebo}} = \left[ \begin{matrix} 0.93 & 0.005 & 0.065 \\ 0.123 & 0.872 & 0.005 \\ 0.026 & 0.031 & 0.942 \end{matrix} \right] \quad \quad \mathbf{A}_{\text{medication}} = \left[ \begin{matrix} 1 & 0 & 0 \\ 0.238 & 0.761 & 0.001 \\ 0.073 & 0.331 & 0.596 \end{matrix} \right]. \end{align*} $$

These estimates are close to those of the MAR model, indicating little initial difference between the conditions, but effectiveness of the drugs reflected in the transition probabilities, which are higher toward the less severe states in the medication compared to the placebo condition.

Results of the state-dependent models for missingness are provided in Table 9. For all three states, the confidence interval for the effect of main excludes 0, indicating a significantly lower proportion of missing ratings at the main measurement occasions. In state 1 and 3, the confidence interval for the effect of week also excludes 0, indicating a higher rate of missing ratings over time, possibly due to attrition. For state 2, the effect of week is not significant. Figure 5 depicts the predicted probability of missing ratings for each state and week. This shows that in state 2, the chance of missing data on the main measurement occasions is small at $p(M_{i,t}|S_{i,t} = 2)=0.003$ , while it is high at $p(M_{i,t}|S_{i,t} = 2) = 0.997$ on the other weeks. In the other states, the probabilities are less extreme, with missing (and non-missing) data occurring on both the main measurement weeks as well as the other weeks. In the final week 6, those in the most severe state 3 are the most likely to have missing data with $p(M_{i,6}|S_{i,6} = 3) = 0.585$ . For those in the least severe state 1, the probability of missingness in week 6 is also substantial at $p(M_{i,6}|S_{i,6} = 1) = 0.272$ .

Table 9

Parameter estimates of the state dependent logistic regression models for missingness, with lower and upper reflecting the lower and upper bounds of the approximate $95\%$ confidence intervals

Figure 5

Predicted probability of missing IMPS Item 79 ratings by week for each state in the three-state MNAR hidden Markov model.

The intercepts for the state-dependent missingness model are also worth considering, especially in interaction with the time-dependent effects. The probability of missingness differs between the states, such that the more extreme states have more missingness. The middle state with medium severe symptoms shows a particular missingness pattern where all the main measurements are almost certainly present whereas the non-main measurements are all missing (see also Figure 6). Both in the least and most severe states, week and main have significant effects. The pattern of correlation between missingness and severity is however complex. In the least severe state, missingness is relatively high at the start and increases only minimally during the study’s 6 weeks. In the most severe state, this is very different: early on the probability of missingness is very low but then steeply increases such that by week 6 the probability of missingness is $p(M_{i,6}|S_{i,6} = 3) = 0.585$ .

Figure 6

Proportions of MAP state assignments over weeks for the medication and placebo groups, according to the MAR and MNAR model.

Disregarding the modeling of missingness, the parameters of the MAR and MNAR model seem reasonably close. This could be an indication that missingness is independent of the hidden states and data are possibly MAR. As discussed previously, the likelihood of the MAR is not directly comparable to that of the MNAR model, as the latter is defined over two variables (the imps79 rating and the binary missing variable), while the former involves just a single variable (imps79). We therefore compare the MNAR model to a constrained version where the parameters of the missingness model are forced to be identical over the states. Unlike the MAR model, this restricted version of the MNAR model accounts for patterns of missingness, allowing these to depend on week and main, but crucially not on the hidden state. A likelihood ratio test indicates that this restricted model fits significantly less well, $\chi ^2(6) = 337.66\ p < .001$ . Hence, there is evidence that the MNAR model is preferable to the MAR model and that missingness is indeed state-dependent.

Whilst the MAR and MNAR model provide roughly equivalent parameters for the severity ratings in the three states, when comparing the MAP state classifications by the Viterbi algorithm (Figure 6), we see that state classifications for the MAR model tend toward the more severe states. According to the MNAR model, during the main measurement occasions, missing values are relatively likely in the least severe state 1. Hence, those with missing values are more likely to be assigned to the least severe state. This is in line with the analysis of Hedeker and Gibbons (Reference Hedeker and Gibbons1997), who found evidence that dropouts in the medication condition showed more improvement in their symptoms before dropping out than those participants who completed the study.

It is worthwhile to note that the MAP states are also determined for time points with missing data, as the transition probabilities make certain states more probable than others, even when there is no direct measurement available. This provides a potentially meaningful basis to impute missing values with e.g., the state-conditional mean. Another option is to impute with an expected rating computed as a weighted sum of all the state-conditional means, weighted by the posterior probability of the states. As imputation is not the focus of this study, we leave the usefulness of such approaches to be investigated in future work.