{'fraction_non_alphanumeric': 0.03916165059105084, 'fraction_numerical': 0.01283631880484444, 'mean_word_length': 4.293474088291747, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A discrete time event-history approach to informative drop-out in multivariate latent Markov models with covariates

Latent Markov (LM) models represent an important tool of analysis of longitudinal data when response variables are affected by time-varying unobserved heterogeneity, which is accounted for by a hidden Markov chain. In order to avoid bias when using a model of this type in the presence of informative drop-out, we propose an event-history (EH) extension of the LM approach that may be used with multivariate longitudinal data, in which one or more outcomes of a different nature are observed at each time occasion. The EH component of the resulting model is referred to the interval-censored drop-out, and bias in LM modeling is avoided by correlated random effects, included in the different model components, which follow a common Markov chain. In order to perform maximum likelihood estimation of the proposed model by the Expectation-Maximization algorithm, we extend the usual backward-forward recursions of Baum and Welch. The algorithm has the same complexity of the one adopted in cases of non-informative drop-out. Standard errors for the parameter estimates are derived by using the Oakes' identity. We illustrate the proposed approach through an application based on data coming from a medical study about primary biliary cirrhosis in which there are two outcomes of interest, the first of which is continuous and the second is binary.


Introduction
In longitudinal studies, subjects may be lost to follow-up due to events, like death, which are associated with the outcome of interest. In these cases, an informative drop-out arises that must be properly modeled in order to avoid biased estimates. From the reverse perspective, the time trend of a longitudinal measurement may predict the risk of an event (e.g., a steadily decreasing CD4 count is predictive of adverse events in HIV patients); see for instance Follmann and Wu (1995) for a general account of related longitudinal and survival processes.
A common approach to deal with informative drop-out is via shared-parameter models (e.g., Wu and Carrol, 1988;Follmann and Wu, 1995), where both longitudinal and survival mechanisms are assumed to share a latent Gaussian variable. Semiparametric shared-parameter models have been proposed by Tsonaka et al. (2009). A discrete random effect, along the lines of this work, is adopted by Roy (2003) to deal with an ordinal latent class model. Another approach to modeling informative drop-out is that of Wulfsohn and Tsiatis (1997) and Rizopoulos (2010), where the risk of an event at time t is influenced by the expected value of the longitudinal response at the same time. The resulting Joint Model (JM) includes both the longitudinal fixed and the random effects in the computation of the hazard function. There are very few generalizations of JMs to the case of discrete longitudinal outcomes. Notable exceptions are that of Rizopoulos and Ghosh (2011), who propose generalized linear models for the longitudinal outcome in a Bayesian framework, and Viviani et al. (2013), who deal with a similar model in a classical maximum likelihood framework.
A limitation of shared-parameter models and JMs is that latent variables, in the form of subjectspecific parameters, are time constant. The effect of unobserved heterogeneity, therefore, cannot evolve or must evolve over time along a pre-specified parametric form (e.g., linear through inclusion of a random slope). Further, the use of continuous latent variables may hinder possibilities of clustering subjects with respect to common unobserved heterogeneity and risk of event/drop-out. Finally, when the outcome of interest is categorical, it may be viewed as measuring, with error, a latent discrete rather than continuous variable. In order to overcome these limitations, latent Markov (LM) models represent a very flexible and convenient way of modeling categorical outcomes which are repeatedly measured over time; see Bartolucci et al. (2012) for an overview. The basic assumption of these models is that the response variables, which are longitudinally observed, are conditionally independent given a hidden first-order Markov chain which accounts for the unobserved heterogeneity.
Despite the relevance of LM models, there are very few extensions of these models to deal with drop-out. A model that is somehow in spirit to our approach is that of Albert (2000), which jointly models the outcome and missing mechanism. The latter is assumed to follow a manifest first-order Markov chain, and the two processes are linked since the outcome is used to model the missingness indicators. In a Bayesian framework, Spagnoli et al. (2011) propose a simple three-state latent chain, in which one state is actually not latent and corresponds to drop-out. Finally, in Maruotti (2013) the time to drop-out is used in a model for the initial and transition probabilities of the hidden Markov chain.
In this paper, we propose a different approach with respect to the ones mentioned above and that, at least to our knowledge, has not been previously considered in the literature. In the proposed approach, the manifest distribution is jointly referred to the longitudinal and drop-out processes; the corresponding time-varying unobserved heterogeneity structure evolves according to the same initial and transition distributions. Our approach falls into the class of non-ignorable random-coefficientbased drop-out models as defined in Little (1995).
In the proposed approach, the longitudinal outcomes are modeled through generalized linear mixed effects models (e.g., McCulloch and Searle, 2004;Fitzmaurice et al., 2004) and a discrete event-history (EH) model (see Steele, 2011, and references therein) is used for the drop-out process. It is important to stress that this is an approach to multivariate longitudinal data, in the sense that even more outcomes and of a different nature can be observed at each time occasions. These outcomes are assumed to be conditionally independent given the random effects. A flexible dependence structure is obtained as the random effects are distributed according to a single first-order homogeneous latent Markov chain with a finite number of states. Subjects in the same latent state share class-specific intercepts for the longitudinal models, and a class-specific intercept for the EH model, and can also share common regression coefficients for the covariates. The resulting estimates are easily interpretable, and the model is reasonable as it is natural to expect that longitudinal outcomes and drop-out share the same sources of (time-varying) unobserved heterogeneity, which can have different effects on each of them.
For the proposed model we perform maximum likelihood estimation by an Expectation-Maximization (EM) algorithm (Baum et al., 1970;Dempster et al., 1977). This requires an extension of the forwardbackward recursions (Baum et al., 1970;Welch, 2003;Bartolucci et al., 2012) to account for informative drop-out. We also pay attention to the computation of the standard errors for the parameter estimates by employing a result due to Oakes (1999); see also Bartolucci and Farcomeni (2013). More-over, for model selection, and in particular for the choice of the number of latent states, we suggest the use of the Bayesian Information Criterion (BIC, Schwarz, 1978).
The remainder of the paper is organized as follows. In Section 2 we illustrate the proposed class of models. In Section 3 we illustrate maximum likelihood inference for these models. Finally, the approach is illustrated in Section 4 through an application based on primary biliary cirrhosis data in which there are two outcomes of interest for each time occasions, one of which is continuous and the other is binary. We provide some concluding remarks in Section 5.

A class of latent Markov models with informative drop-out
We consider a longitudinal study on a sample of n subjects, or more generally sample units, in which s follow-up time occasions are scheduled for each subject. We assume that s is known in advance and equal for all subjects. A generalization to a subject-specific number of occasions is straightforward, while if s is not known in advance all proposed inference can be thought of as being conditional on the maximum number of follow-up times which is observed.
For every subject i and time occasion t, with i = 1, . . . , n and t = 1, . . . , T , we observe r response variables, denoted by Y hit , h = 1, . . . , r; we also denote by T i the last time occasion of observation for subject i, so that drop-out occurs before occasion T i + 1. As will be clear in the following, also T i is a random variable, the distribution of which depends on observable covariates and latent variables for the unobserved heterogeneity. A realization of the h-th response variable is denoted by y hit and, accordingly, a realization of T i is denoted by t i . The observed outcomes for the same subject i and time occasion t are collected in the column vector y it = (y 1it , . . . , y rit ) ′ . Also note that if there is drop-out then t i < s, whereas t i = s indicates that a complete record of outcomes is observed for the i-th subject.
Let D it , i = 1, . . . , n, t = 1, . . . , s, denote a binary random variable equal 1 if subject i drops out from the study after occasion t and before occasion t + 1, that is, T i = t, and to 0 otherwise. The basic assumption of the proposed model is that, given a discrete latent variable U it with k support points and the available covariates, the response variables Y 1it , . . . , Y rit are conditionally independent and they are also independent of D it . In particular, we denote by x hit the column vector of covariates affecting Y hit and by µ hit (u) the expected value of this response variable given these covariates and U it = u. Note that the dependence on the covariates of this mean is not explicitly indicated since these covariates are considered as fixed and given; this convention will be used for the notation throughout the paper. Similarly, we denote by z it the column vector of covariates affecting D it and by p it (u) the probability that D it = 1 given these covariates and U it = u. Then, for u = 1, . . . , k, we assume where g h (·) and g(·) are appropriate link functions and every Y hit is assumed to have conditional distribution belonging to the regular exponential family (McCullagh and Nelder, 1989). The reason why p is (u) = 1 is that the follow-up surely stops after occasion s for all subjects.
In the dataset used to illustrate the proposed approach (see Section 4) there are r = 2 response variables. The first of these variables is continuous and the second is binary. Hence, for this application we choose g 1 (·) as corresponding to the identity link and g 2 (·) as corresponding to the logit link; for the distribution of the first variable we also have a dispersion parameter indicated, in general, by σ 2 h . The model for each longitudinal outcome is a classical generalized linear mixed effects model, while the time to drop-out follows a geometric distribution as in classical discrete time EH models (e.g., Steele, 2011). We regard drop-out as a trial within each time interval. It is shown in Steele et al. (1996) that the resulting likelihood is that of a Bernoulli; consequently, we have where, since p is (u) = 1, the probability p it i (u i ) disappears from the above expression when there is no drop-out until the end of the study (t i = s). This recovers the truncated geometric distribution given that s is finite. In order to derive the results in Section 3 it is also important to note that and that Note that the only information that is used in the model is the interval censored event time (i.e., that drop-out occurs between t i and t i+1 ). We further note that when continuous-time durations are grouped into discrete intervals, a continuous-time hazard model (with constant hazard within each interval) would lead to the complementary log-log model with g[p it (u)] = log{− log[1 − p it (u)]}; see Kalbfleisch and Prentice (2002). In our application we use a logit link, which we find more convenient.
Regarding the distribution of the latent variables, we assume that, for i = 1, . . . , n, the sequence U i1 , . . . , U is follows a Markov chain with initial probabilities which are collected in the column vector π, and time-homogeneous transition probabilities which are collected in the transition matrix Π. The parameters of the model that must be estimated are then π, Π, α, β, γ, δ, where α and γ are column vectors with elements α hu , h = 1, . . . , r, u = 1, . . . , k, and γ u , u = 1, . . . , k, respectively. Similarly, β is the vector obtained by casting the vectors β h , h = 1, . . . , r.
The degree of dependence between the longitudinal and the survival processes is measured by the total variation of the support points of the EH model. Consequently, when k = 1 then drop-out is non-informative. A classical LM model with k > 1 states, where drop-out is non-informative, may obviously also be obtained by fixing γ u = γ, u = 1, . . . , k. There is sensitivity to drop-out as soon as γ is not constant, hence we can use a formal test based on a likelihood ratio statistic for the null hypothesis H 0 : γ 1 = · · · = γ k .
Note that the link between the random effects γ u and α hu is based on the assumption that they follow the same latent Markov rule. This approach is similar in spirit to situations in which a copula is used to model the dependence of random effects of the two processes, like in Rizopoulos et al. (2008).

Likelihood inference
We begin considering the observed likelihood where θ is a short-hand notation for all the model parameters and f (y i1 , . . . , y it i , t i ) is the density or probability of the observed outcomes given the covariates x it and z it until occasion t i ; these covariates are not explicitly indicated as they are fixed and given. This component of the observed likelihood can be expressed as where f (y hit |u) refers to the conditional density or probability of Y hit evaluated at y hit , given U t = u the corresponding covariates, and , for t = 1, . . . , s − 1 with p is (u) = 1.
Expression (1) can be efficiently computed by an extension of the forward recursion, which is well known in the hidden Markov literature (Baum et al., 1970;Zucchini and MacDonald, 2009;Bartolucci et al., 2012). First of all, ruling out the trivial case in which s = 1, we consider the following density or probability for t = 1, . . . , t i and u = 1, . . . , k: Then, for t = 1 we have that whereas, provided that t i > 1, for t = 2, . . . , t i we have where v = 1, . . . , k. At the end of the recursion (t = t i ), we have that In order to maximize the likelihood L(θ), we implement a version of the EM algorithm (Dempster et al., 1977), which is based on the complete data likelihood. Let w it (u) denote a dummy variable equal to 1 if the i-th subject is in latent state u at the t-th occasion and let z it (u, v) = w i,t−1 (u)w it (v) be a dummy variable equal to 1 if there is a transition from latent state u to latent state v at occasion t.
The logarithm of this likelihood has the following expression: where the second summand involving z it (u, v) disappears if t i = 1.
With reference to our application based on r = 2 response variables, the first of which has normal conditional distribution parametrized by an identity link function and the second of which has Bernoulli distribution parametrized by a logit link function, we have that where φ(·) denotes the standard normal density function.
The EM algorithm alternates two steps until convergence: first, the conditional expected value of the complete data log-likelihood is obtained (E step). The resulting expression is then maximized with respect to θ (M step).
The EM is guaranteed to converge to a local optimum of the observed likelihood. In order to increase the chances of reaching the global maximum, we use a multistart strategy. The first initialization is based on estimating β and δ using separate generalized linear models. The parameters α and γ are then centered on the intercepts obtained by the previous generalized linear models. Finally, Π is initialized so that all off-diagonal elements π uv , with v = u, are equal to 1/[k(k − 1)]. The other starting solutions are obtained by randomly perturbing the parameter estimates obtained at convergence from the deterministic one.
Finally, we use the value of the likelihood at convergence to compute the index on which the Bayesian Information Criterion (BIC) for model choice is based (Schwarz, 1978).

E step
At the E-step, the conditional expected value of (2) is simply computed by a plug-in of the expected values of w it (u) and z it (u, v) given the observed data and the current value of the parameters. These expected values are denoted by w it (u) and z iu (u, v), and can be computed by means of an appropriate backward recursion adapted from the hidden Markov literature as we illustrate below. First of all, we consider the probabilities: We have that b it i (u) = 1 for u = 1, . . . , k, whereas, for t = 1, . . . , t i − 1 and u = 1, . . . , k, we have where as before p is (u) = 1. At the end of the recursion, we obtain the following expected values: , t = 1, . . . , t i , u = 1, . . . , k; provided that t i ≥ 2, we also obtain, for u, v = 1, . . . , k, the expected values

M step
At the M-step the conditional expected value of (2) is maximized by separately maximizing its addends.
It is straightforward to check that explicit solutions are available for π u ; in particular we have In a similar way, for the transition probabilities we have the following solution: where the outer sum at numerator and denominator is extended to all i such that t i ≥ 2.
For what concerns the other parameters, we can use r + 1 separate Newton-Raphson algorithms, similar to that used for standard generalized linear models. The first r algorithms are used to update the parameters α hu and β h by maximizing the following expression: for h = 1, . . . , r. For instance, in the case of a binary outcome with logit link, this corresponds to the maximization of whereas explicit solutions are available with Gaussian outcomes also for the dispersion parameter σ 2 h . The last Newton-Raphson is used to maximize, with respect to γ and δ, the expression

Computation of standard errors
In order to estimate the standard errors for the parameter estimates, we use an approach based on Oakes (1999)'s equality, simplified by the fact that the score vector for the expected complete loglikelihood is often available in a closed form, given that it corresponds to the score of r + 1 generalized linear models; see also Bartolucci and Farcomeni (2013). For instance, for binary outcomes with logit link, the first derivative of the expected complete log-likelihood with respect to ξ hu is easily derived as: whereas the first derivative with respect to β h can be derived as Moreover, the first derivative of the expected complete log-likelihood with respect to γ u can be derived as: and finally the derivative with respect to δ is These derivatives are the same that are used in the Newton-Raphson algorithm to implement the M-step as described above.
Once we have the score vector, it can be noted as in Oakes (1999) that the observed information is equal to the Jacobian of the score vector with respect to the parameters, plus the Jacobian of the score vector with respect to w it (u), seen as a function of the parameters. For convenience, the second Jacobian is computed numerically.

Application to primary biliary cirrhosis data
We illustrate the proposed approach by an application based on data coming from a randomized study for treatment of primary biliary cirrhosis. These data were previously analyzed by Murtaugh et al. (1994) and Rizopolous et al. (2010) from a slightly different perspective than the present one. In this study, n = 312 patients were randomized to a placebo or a treatment based on D-penicillamine. We are interested in evaluating the effect of treatment on a continuous outcome (serum Bilirubin in mg/dl, Y 1 ) and a binary outcome (presence of edema, Y 2 ), after adjusting for certain covariates (drug, age, gender, albumin in gm/dl, alkaline phosphatase in U/liter, and SGOT in U/ml at the first visit) and drop-out. Continuous covariates are zero centered, as the sample mean has been subtracted from each of these covariates. In this application it is reasonable to expect drop-out to be informative as it may be due to death related to the illness or to a transplant.
The maximum number of follow-up time occasions is s = 16 for these data and only 1% of patients have a complete record; the median time to drop-out is 5. Note that a high serum Bilirubin is very likely to speed the cirrhosis up, hence making drop-out due to death more likely; this enforces our idea that the drop-out cannot be considered as non-informative. In this regard, Table 1 reports the proportion of subjects having a certain number of observations, that is, n i=1 I(T i = t)/n, t = 1, . . . , T , the corresponding Kaplan and Meier (1958)'s estimates, and the mean of the two outcomes based on the number of survivors.
[ Table 1 about here.] We begin the analysis of these data by selecting the number of latent states of the hidden Markov chain. In Table 2 we report the maximum log-likelihood, number of parameters, and BIC (Schwarz, 1978) for increasing values of k from 1 to 4.
[ On the basis of the results in Table 2, we select k = 2, as the corresponding model has the lowest BIC. In fact, for this model we have maximum log-likelihood of -7265.3 with 34 parameters and then BIC=14725.9. Sensitivity to drop-out can be evaluated by fitting the same model with k = 2, but with constant latent intercepts (i.e., H 0 : γ 1 = γ 2 ). The latter model has a maximum log-likelihood of -7334.5 with 33 parameters and then BIC=14858.5, which is much higher than the previous value.
Then, we do not consider hypothesis H 0 to be plausible at least when k = 2, as the likelihood ratio test statistic for this hypothesis, which is equal to 138.4, confirms. The conclusion is that there is evidence of informative drop-out with these data.
In Table 3 we report the parameter estimates obtained with our model, and the corresponding LM model that assumes the drop-out to be non-informative (including H 0 ). In this table, the intercepts are obtained by averaging the support points.
[ Table 3 about here.] From the left panel of Table 3 (results under the proposed model) we can draw the following conclusions. First of all, the treatment does not seem to be effective on time to drop-out or on level of serum Bilirubin, but at least it significantly decreases the probability of edema, with an OR of 0.57 after adjusting for other factors. Note that, instead, Rizopolous et al. (2010) concluded that the treatment is effective on serum Bilirubin. In our interpretation, this difference is due to the different way of modeling unobserved heterogeneity, that in the approach of Rizopolous et al. (2010) is assumed to increase linearly in time. As far as the other predictors are concerned, it can be seen that with age the level of serum Bilirubin decreases and the probability of edema increases, that females are at higher risk of drop-out and edema but have lower levels of Bilirubin, and that albumin is protective.
A comparison between the left panel and the right panel of Table 3 allows us to study the sensitivity to drop-out. Assuming that the drop-out is non-informative (under H 0 ) leads to certain differences in the parameter estimates and standard errors, especially for the time-event history part for which we may observe that several parameter estimates change sign under H 0 with respect to the initial model specification.
We finally consider the estimated latent distribution parameters. In Table 4 we report the differences between the latent intercepts and their averages (which are reported in Table 3), the initial parameter vector, and the transition matrix.
[ We observe that the two groups are very well separated. Three out of four patients start in group 1, which has a slightly lower serum Bilirubin, propensity to edema and probability of drop-out. This group is highly persistent, but still it can be seen from the transition matrix that the proportion of subjects in this group steadily decreases over time. At t = 5, only 60% are in this group, and at t = 15 this proportion has decreased to 45%. As a consequence, it can be said that the global health status of patients has worsened over time more than as predicted by the trend estimated by fixed effects (Table   3). In particular, the second group is made instead of patients with a sensibly higher Bilirubin, much higher propensity of edema and slightly higher risk of drop-out.

Conclusions
In this paper, we propose an event-history approach (Steele, 2011) to account for informative dropout in LM models for multivariate longitudinal outcomes. The main features of these models are that they assume conditional independence among response variables given a hidden Markov chain. In order to deal with the models formulated in this way, we extend the usual forward-backward recursions (Baum et al., 1970;Dempster et al., 1977) to informative drop-out and consequently propose an extended version of the Expectation-Maximization algorithm (Dempster et al., 1977) for maximum likelihood estimation of LM models.
One of the main features of an LM model with covariates is that it allows for a form of unobserved heterogeneity that is not restricted to be time constant. This form of heterogeneity is typically explained on the basis of missing covariates, and its evolution is accounted for by the inclusion of random intercepts in the model. These intercepts are time and individual specific and follow a Markov chain.
Assuming time-constant unobserved heterogeneity may be restrictive, especially when individuals are followed up for many occasions and/or for a long period of time. The advantages of the proposed type of modeling are shown by an application based on a study about primary biliary cirrhosis in which there are two outcomes of interest, the first of which is continuous and the second is binary.
An important assumption of the adopted LM model is that the response variables referred to the same time occasion are conditionally independent given the corresponding latent variable and covariates; this assumption is usually known as local independence. This assumption simplifies model estimation substantially and may be explained by considering that all factors affecting the responses, that is, observed and unobserved covariates, are properly accounted for. In any case, this form of local independence can be relaxed by using a multivariate link function, as the one adopted in the LM framework proposed by Bartolucci and Farcomeni (2009). Table 1: Observed mean of serum Bilirubin and proportion of subjects with edema by time for the primary biliary cirrhosis data, together with the proportion of survivors ( n i=1 I(T i = t)/n) and corresponding Kaplan and Meier (1958)'s estimates (KM t Table 4: Differences between the latent intercepts and their average, initial parameter vector and hidden transition matrix for the primary biliary cirrhosis data obtained with the proposed LM model. The number of latent states is fixed at k = 2.