Partially hidden multi-state modelling of a prolonged disease state defined by a composite outcome

For rheumatic diseases, Minimal Disease Activity (MDA) is usually defined as a composite outcome which is a function of several individual outcomes describing symptoms or quality of life. There is ever increasing interest in MDA but relatively little has been done to characterise the pattern of MDA over time. Motivated by the aim of improving the modelling of MDA in psoriatic arthritis, the use of a two-state model to estimate characteristics of the MDA process is illustrated when there is particular interest in prolonged periods of MDA. Because not all outcomes necessary to define MDA are measured at all clinic visits, a partially hidden multi-state model with latent states is used. The defining outcomes are modelled as conditionally independent given these latent states, enabling information from all visits, even those with missing data on some variables, to be used. Data from the Toronto Psoriatic Arthritis Clinic are analysed to demonstrate improvements in accuracy and precision from the inclusion of data from visits with incomplete information on MDA. An additional benefit of this model is that it can be extended to incorporate explanatory variables, which allows process characteristics to be compared between groups. In the example, the effect of explanatory variables, modelled through the use of relative risks, is also summarised in a potentially more clinically meaningful manner by comparing times in states, and probabilities of visiting states, between patient groups.


Introduction
For studies in rheumatic diseases, and in other medical contexts, the outcome variable of interest is often composite. Such an outcome will be defined based on the observed values of a set of separate variables that all reflect some aspect of a patient's disease activity. Sometimes the composite outcome is a mathematical function of the values of the constituent variables and sometimes it may be a categorical variable representing disease states defined in terms of the constituent variables. In this paper, we focus on the latter situation.
It may also be the case that clinical interest focuses on a patient being in a disease state for a prolonged period of time. For example, the concept of minimal disease activity (MDA) in rheumatic disease was conceptually defined as "that state of disease activity deemed a useful target of treatment by both the patient and physician, given current treatment possibilities and limitations" by The Outcome Measures in Rheumatology Clinical Trials 6 Conference (Wells et al. 2005). This reflects the fact that the complete absence of disease is not a realistic goal for many patients. For psoriatic arthritis (PsA), an inflammatory arthritis associated with the skin disease psoriasis, this has been operationally defined in terms of 7 criteria related to physician, patient and laboratory measures of disease activity, that is, disease symptoms that are potentially reversible through treatment or other factors. However, short term MDA is of little clinical interest as it is MDA of extended duration, typically one year or more, that has been linked with reducing the risk of permanent joint damage, a major aspect of disease progression in PsA (Coates et al. 2010a).
There are challenges to the analysis of events that are defined by prolonged observation of a condition (Farewell and Su 2011) and relatively simple approaches are often adopted in practice. For example, Coates et al. (2010a) divided a longitudinal cohort of patients into two groups, those who achieved the criteria for MDA at consecutive visits for a minimum of 12 months and those who did not over their periods of followup. These two groups were then compared in various ways in terms of explanatory variables. This approach does not appear to take full advantage of the longitudinal follow-up of the patients or reflect the intermittent observation patterns of the cohort. Along with the need for more comprehensive longitudinal modelling reflecting intermittent observation of patients, typically at clinic visits, a sizeable number of clinic visits may not provide information on a sufficient number of MDA criteria to unambiguously determine whether a patient is in the MDA state. In this paper, we examine how these challenges may be met when adopting a simple two-state model for the presence and absence of MDA in PsA. The primary aim is to provide a means to characterise the MDA process in PsA. An additional benefit is that this model can be extended to incorporate explanatory variables, which allows process characteristics to be compared between groups. The effect of explanatory variables, modelled through the use of relative risks, can also be summarised in a potentially more clinically mean-ingful manner by comparing times spent in states, and probabilities of visiting states, between patient groups.
As in Sweeting et al. (2010), features of the observation process suggest the use of a partially hidden multi-state model. Aalen (2010) commented on this previous work that the model formulation was needed to address the nature of the available data. He commented on the use of Markov formulations as "a simple way of introducing dynamics into the system" and that while "the Markov assumption is often criticized as being too strong … a simple Markov assumption will, for many purposes be good enough". In addition, our incorporation of available data from visits when MDA can not be unambiguously determined, is also similar to the use of an auxiliary variable in Sweeting et al. (2010) to address the problem of informative observation, another consideration highlighted in Aalen (2010). Thus the modelling approach discussed in this paper is, we hope, taking account of the issues raised in Aalen (2010) and it is a great pleasure to contribute the paper to this special issue of the journal prepared in honour of Odd Aalen's long and distinguished research career.

The clinical example
Our motivating example is based on data from 7024 clinical visits from 856 patients seen at the University of Toronto PsA Clinic since 2003. Patients were evaluated using a standard protocol every 6-12 months. Patients were followed up for a median time of 3 years (maximum 10 years), with a median 6 visits (maximum 27). This intermittent observation pattern needs to be reflected in analyses, as discussed earlier, but it is important to note that visits occurring within 3 months of a regularly scheduled visit, perhaps to address clinical needs identified at the previous visit, are not included in the database. Clinical assessments included the number of (out of 68) tender joints and the count of (out of 66 excluding hips) swollen joints, a measure of enthesitis reflecting the number of inflamed locations where tendons or ligaments insert into bones, and a dactylitis score reflecting the extent to which entire digits are inflamed, a characteristic symptom of PsA. Skin assessment included both the body surface area (BSA) and the Psoriasis Area and Severity Index (PASI) which has a range 0-72. A clinically count of permanently damaged joints was also recorded at each visit. A physician global assessment on a 10 cm scale was to be completed at every visit and patients completed self-reported questionnaires including the Health Assessment Questionnaire (HAQ), which has a range of 0-3, and patient global assessments, on a 10 cm scale, usually at every other visit.
The criteria for the definition of MDA used by Coates et al. (2010b) were 5 or more out of 7 of the following: 1. Tender joint count (TJC) ≤1 2. Swollen joint counts (SJC) ≤ 1 3. PASI score ≤1 or BSA ≤ 3% 4. Patient pain visual analog score (PTPPAINV) ≤ 1.5 cm 5. Patient global disease activity visual analogue score (PTPSA) ≤ 2 cm 6. HAQ score ≤ 0.5  All but 8 clinic visits provided information on at least one of the MDA criteria, and the numbers with 0 to 7 criteria missing were 1367,2807,1449,924,357,96,16 and 8 respectively. The HAQ and patient global scores were missing approximately 50% of the time as per their scheduled administration at every other visit and the rest were missing at 10% to 20% of visits. The MDA status could be determined at 63% of the visits comprising 38% when at least 5 out of 7 criteria were observed and satisfied and 25% when at least 3 out of 7 were observed and not satisfied. There were 1390 and 195 occasions when MDA was observed followed by no MDA and MDA, respectively, at the next visit, and there were 144, and 825 occasions when no MDA  Table 1 presents the number of positive criteria observed by MDA status and the number of criteria observed. For the purposes of regression analyses, in which the effect of baseline explanatory variables on MDA is subsequently examined, if no treatment information is recorded for a patient then it has been assumed that neither disease modifying anti-rheumatic drugs (DMARDs) or biologic agents were given. Patients with missing information on any baseline explanatory variable were excluded. The highest fraction of missingness, 9%, was seen with the binary indicator for the involvement of axial joints. The other indicator variables used in our example analyses were for polyarthritis (the involvement of five or more joints), sex, an elevated sedimentation rate (ESR) and previous damaged joints. Age and disease duration prior to clinic entry were also included as continuous variables. For each patient, let S j represent their MDA status at clinic visit j, and let y jk represent the observed value of the kth variable used to define MDA status at clinic visit j. Because two variables may be used to determine the third MDA criterion described in Sect. 2, there are 8 defining variables in total. Additionally, let x j represent a row vector of explanatory variables associated with the patient at visit j. The transition rates can then be specified as

The model
and where λ (No→MDA,0) and λ (MDA→No,0) are baseline intensities, β and γ are column vectors of regression coefficients associated with the explanatory variables in the two models. Note that the model is specified in continuous time although observation is intermittent at clinic visits.
Equations (1) and (2) reflect the necessary simplifying assumptions for a timehomogeneous Markov multi-state model with constant transition intensities, although this may be relaxed. The time-homogeneity assumption can be relaxed easily with available software through the use of piece-wise constant transition intensities. This might be needed, for example, if the introduction of new treatments resulted in variation of MDA frequency with calendar time. Departures from the Markov assumption would introduce more complications since state entry times are unknown, though some non-Markov models can be fitted to data of this kind in available software using phase-type sojourn distributions (Titman and Sharples 2010). This is expanded upon in Sect. 5. The Markov assumption is particularly useful for fitting a partially hidden multi-state model and for calculating summary characteristics of the MDA process. It is not expected that the Markov assumption which introduces the dependence of the future on the past through conditioning on current state would be an undue simplification in the context of PsA. With respect to the use of this model, in comparison to the very simple models used previously (Coates et al. 2010b), we would regard the multi-state model structure to be a critical assumption to adequately model the MDA process, but the specific Markov assumption to be a simplifying assumption regarding a secondary aspect of the model, following the approach to assumptions outlined by Cox and Snell (1981).
Because MDA status can not be determined at all visits, it is convenient to regard this model as a partially hidden multi-state model. At some visits the MDA status is known but at others it is unknown or hidden. This essentially extends the usual multistate modelling approach to allow information from the y jk variables to provide extra information on MDA status at visits when the binary classification of MDA based on the y jk variables cannot be unambiguously determined. This could be done by incorporating modelling of the conditional distributions of the binary criteria derived from the y jk variables, given MDA status, but we will focus only on the somewhat more general approach of directly modelling the conditional distributions of the 8 y jk variables that determine these binary criteria. A comparison of these two approaches can be found in Jackson et al. (2016) where there is some evidence that modelling the y jk variables can provide greater precision for estimation of parameters in the multi-state model, as might be expected.
It is assumed that given the (observed or latent) MDA status, the distributions of the y jk variables are independent from each other. In other words, we assume that the marginal distributions of y jk variables help to discriminate the MDA status, but the associations between y jk variables will not provide additional information for this discrimination. Without this conditioning, an independence assumption would be unreasonable but it is less problematic given the conditioning, although it is unlikely to be exactly true.
In terms of missing data for y jk , we assume that the unobserved y jk variables are missing at random given the observed y jk (k = k ) values at visit j. Therefore, we don't model missing indicators of y jk and relate them to the latent MDA status. If missingness depends on the unobserved y jk values after conditioning on the observed data, then the missing indicators should also inform the latent MDA status and need to be modelled. This will correspond to a latent ignorability assumption discussed in Harel and Schafer (2009). Because a substantial amount of partially missing data in the PsA clinic are due to different schedules for data collection, e.g., HAQ is only measured every other visit, we reckon that the missing at random assumption is reasonably plausible in this context.
To specify the probabilities Pr(y jk |S j = r ), the patient pain and global activity scores are rounded to integers and assumed to arise from Binomial(10, p kr ) distributions while the remaining variables, which are all integers if HAQ and PASI are multiplied by 100, are specified to arise from negative binomial distributions, N eg Bin(n k , p kr ). Heuristically, information on the set of all n k and p kr parameters will arise primarily from visits when S j is observed, while for latent or hidden values of S j , the subset of y jk values observed will provide information on the possible MDA status. Low values of y jk variables are more likely to be associated with an underlying MDA state.

Estimation
The proposed model, with its Markov assumption, can be fitted by full maximum likelihood. Introducing an additional subscript i for patients, let y i j represent the vector of MDA defining variables observed at visit j from patient i, where j = 1, . . . , n i and i = 1, . . . , m.
Let Pr(S i j | S i, j−1 , q, x i j ) be the transition probability from state S i, j−1 to state S i j over the time interval separating visits j − 1 and j, given explanatory variables x i j , where q = (λ (No→MDA,0) , λ (MDA→No,0) , β, γ ) represents both the transition rates governing the hidden Markov chain and the effects of explanatory variables on these. Let Pr(S i1 | f) be the distribution of (potentially unknown) MDA states at the initial visit, with vector of probabilities f. Finally, let f (y i j | S i j , α) be the conditional distribution of y i j given the states S i j = 0 ("no MDA") and S i j = 1 ("MDA"), with the parameter vector α. Then, assuming that the y i j are conditionally independent given S i j , the full likelihood can be represented as where the parameters are θ = (q, f, α), {S} is the set of all observed MDA states, and the product over visits j is summed over all possible latent state pathways {S i } for each patient i (Satten and Longini 1996). Note that the "data" in this model implicitly includes the observations of MDA status S i j at times j when this is known, which constrains the set of latent state pathways to be summed over. Satten and Longini (1996) showed further that the likelihood contribution from a patient i in this model can be expressed as a product of n i K × K matrices, where K is the number of states in the Markov model structure (K = 2 in our example), which facilitates computation. Our model generalises the model in Satten and Longini (1996) to composite outcomes given the hidden state, and a combination of observed and hidden states S i j .
Maximum likelihood estimation is implemented in the msm R package (Jackson 2011; R Development Core Team 2010) for continuous-time Markov and hidden Markov modelling. The package allows general state-transition structures with transition intensities depending on explanatory variables. There can be any number of outcomes linked to a hidden state, with a variety of distributional assumptions possible. The implementation is based on derivatives of the log-likelihood (Lystig and Hughes 2002) and uses the R optim BFGS method.

Complete case analysis
If the possible information from y i j values observed at visits when MDA status can not be determined is ignored, then the two-state model in Fig. 1 can simply be fitted using the standard likelihood for a continuous-time Markov model for panel data (Kalbfleisch and Lawless 1985), where Pr(S i j |S i, j−1 , q, x i j ) is the transition probability from observed state S i, j−1 to state S i j over the time interval separating visits j −1 and j given x i j , and the likelihood contribution for each person is conditioned on their initial observation S i1 . This can again be implemented in the msm package. Note that in this model, the observations include only the 63% of patient visits at which MDA can be determined. This analysis will be termed a complete case analysis. The comparison of this with the analysis based on the partially hidden multi-state model analysis will provide some indication of whether the latter can provide any gains in precision or any bias reduction relative to the former.

Analyses related to sustained MDA
As well as parameter estimation of a multi-state model, there is often interest in summary measures related to state occupancy. Estimation of quantities such as the expected duration of time in a state and total time in a state or the number of times in a state over a fixed time period can be derived as analytic functions of transition rates from continuous-time Markov chain theory. However, there is clinical interest in prolonged durations of state occupancy for MDA, such as the one year duration to define sustained MDA. For the estimation of these, simple analytic calculations of relevant measures are not possible.
Therefore, in order to provide information on sustained MDA, estimated expectations can be calculated by simulation of state histories over a period of time, say 10 years, for 100,000 patients under the fitted multi-state model. Specifically, the history of a patient with explanatory variable vector x j is simulated as a series of periods alternately spent in No MDA and MDA, starting with No MDA, and each with duration simulated from an exponential distribution with rate λ No→MDA (x j ) or λ MDA→No (x j ) respectively, until a censoring point of 10 years. Then the probability of visiting sus- tained MDA, for example, can be estimated as the proportion of simulated patients with sojourns of one year or more in the MDA state starting before 10 years. If the patient is in MDA at 10 years, then they are followed up further until the end of the MDA sojourn, so that this sojourn can be categorised as sustained or not. Furthermore, as done in Aalen et al. (1997) for simpler functions of parameters from a multi-state model, confidence intervals and standard errors for these estimates can be determined by simulating from the distribution of parameter values, say 1000 times, and repeating the simulation of state histories for 100,000 patients for each of these 1000 sets of parameter values.

Simple two-state model
Figure 2 (first and second columns) shows the observed and fitted distributions of each MDA-defining criterion conditional on known MDA states along with an indication of the binary cutpoints used to define the binary MDA defining criteria. Not surprisingly, some variables appear to have more potential to discriminate between the two states than others. Also shown in Fig. 2 are vertical bars within each histogram segment which display the maximum likelihood estimation results for these distributions. The negative binomial model can be seen to fit well for TJC, SJC and ENTH_TOT. However for PASI, BSA and HAQ, the shape of the negative binomial does not perfectly represent the spike at zero and the distribution of the non-zero values. Similarly, for PTPAINV and PTPSA, the variance of the observed distributions can be seen to differ slightly from the variance of the fitted binomial. Divergence from the histograms which could result from an inappropriate distributional assumption and/or the fact that the fitted distributions make use of additional data, as shown in columns three of Fig. 2, from patients visits at which MDA can not be determined unequivocally but which generate weighted contributions to the fit of the conditional distributions of interest. To consider the possibility of inappropriate distributional assumptions, an alternative model was explored. For PTPAINV and PTPSA, beta-binomial conditional distributions were fitted, to account for the over/under-dispersion. For the remaining three criteria, to robustify against model misspecification, these were coarsened into two binary criteria, according to the definition of MDA from Coates et al. (2010b), that is, HAQ > 0.5, and a single outcome of PASI ≤ 1 or BSA ≤ 3%, and binary conditional distributions were fitted. The estimate of the expected total years in MDA over 10 years under this model is 4.43, compared to 4.47, as reported subsequently, under the original model. Therefore this key result appears to be robust to specification of the conditional distributions of the outcomes. The first two lines of Table 2 provide estimated mean lengths of one period in the MDA and No MDA states as well as the associated standard errors based simply on the two estimated transition rates from the model without explanatory variables in Fig. 1. Results are provided from both the complete case analysis and the fitted partially hidden multi-state model based on data from all patient visits. It can be seen that there is an increased precision of estimation from the latter but also that the estimated times are substantially less. Thus, there is evidence of potentially notable bias in the complete case analysis. This may arise due to the variation seen between visits in the outcome variables which suggests greater movement between states than would be evident from the complete case analysis with its longer periods between observations. The lower section of Table 2 presents estimation results for the expected total time in MDA, the expected number of MDA periods and the probability of visiting MDA at least once over a 10 year period. As well as providing results from the complete case and the partially hidden multi-state model analyses, estimates are also provided, through simulation as outlined in Sect. 3.4, for only MDA periods which last longer than one year. As expected given the results for the length of times in the states, the more complete use of the available data generates increased estimates for the expected total time in the MDA state, the expected number of periods of MDA and the probability of at least one period of MDA. And, again, as would be expected, these values are all reduced when focus is only on MDA periods of sustained length.
Note that the results in Table 2 are influenced by the 10 year horizon. For example, in the right column, the expected years in MDA of 4.47 is not the product of the average duration in MDA, 3.10, and the expected number of entries, 1.96, because of the 10 year cut-off when some patients would be expected to be in the MDA state.

Explanatory variables
Figure 3 presents estimated transition intensity ratios (HRs), and their 95% confidence intervals, from the multi-state model when baseline explanatory variables, as determined when first entering the cohort, are incorporated into the transition rate functions as specified in Sect. 3.1. The presentation is restricted to the results from two multivariable analyses, using complete cases and using the partially hidden multi-state model, incorporating the variables shown in Fig. 3. It can be seen that the most notable effects are associated with sex and with disease presentation as reflected in polyarthritis and axial joint involvement.
While the parameters represented in Fig. 3 derive from a very convenient relative risk model for the effects of explanatory variables, it is perhaps difficult to communicate the overall clinical implications of these effects. For example, females appear to be less likely to enter MDA and more likely to leave but it is useful to have some indication of how these effects combine to create different patterns of disease. Tables 3 and 4 gives an illustration of how this might be done. For simplicity two single factor partially hidden multi-state model regression analyses, one including a binary indicator of female sex and the other the two binary indicators for polyarthritis and axial joint involvement are examined. These analyses, which do not adjust for other explanatory variables and therefore not conditioned on them, are not directly comparable to that presented in Fig. 3 but comparable calculations could be done for any single factor holding other factors constant using this larger model. Measures of various aspects of state occupancy for these two models are presented in Table 3 and relative measures are presented in Table 4. Some of these measures can be calculated analytically but also given are confidence intervals, all of which are derived from the simulation approach of Sect. 3.4. The more positive prognosis for males in regard to MDA can be clearly seen in the relative measures of Table 4, except for spells in MDA where the number of spells is similar for males and females, being 1.84 for males and 1.86 for females.
For the model including disease pattern variables, the relative size of the two regression coefficients seen in Fig. 3 is reflected in Table 4 by more dramatic effects for axial involvement, or both axial and polyarthritis, relative to neither disease pattern being present. Note that for both models, the effects on the odds ratio scale appear dramatic. However, this derives partially from the high probabilities of at least one MDA spell, for example 0.98 for males and 0.94 for females, so that odds ratios in this probability range can be extreme although the absolute difference in probabilities is small.

Discussion
A partially hidden multi-state model provides a framework for studying intermittently observed composite outcomes such as MDA. Notably, it provides a natural way to incorporate observations from the constituent variables that define a composite outcome from observation times when not all these variables are observed and the composite outcome can not be determined. The analyses presented in this paper for the specific case of MDA in PsA illustrate the potential for this to increase precision and to protect against bias. Coates et al. (2010a) previously examined MDA in psoriatic arthritis but as well as requiring 5 of the 7 criteria to be fulfilled, also required that MDA must be observed at  consecutive visits for a minimum of 12 months in order to focus on sustained MDA. In our example dataset, which updates that of Coates et al. (2010a), and based on complete case data, 229/619 (37%) of patients achieved this and the median duration of such episodes was 42 months (3.5 years), greater than the median of 28 months presented in Coates et al. (2010a) based on earlier data on 344 patients. While there may be other reasons for this difference, the difference is at least partially explained simply on the basis of followup times as the length of MDA episodes will be censored at the last observation time. For these episodes in our data which begin prior to 2008, which is the cutoff for the data of Coates et al. (2010a), the mean duration is 76 months (6.3 years), reflecting the additional followup of the patients considered in Coates et al. (2010a). For MDA episodes in our data beginning after 2007 the mean duration is 27 months (2.3 years). Thus, estimation of the length of MDA episodes in this manner is problematic and the estimated mean durations arising from a two-state model should be preferred as these are valid estimates not influenced by followup times. As a check of the Markov assumption used in the models reported, a semi-Markov model was fitted to the data with fully-observed MDA statuses, using "phase-type" distributions. The two states are divided into two latent "phases", resulting in a fourstate hidden Markov model in Fig. 4, with 6 instead of 2 transition rates to be estimated. Thus the exponentially-distributed sojourn in each state is replaced by a sequence of either one or two sojourns with different transition rates. This allows the transition intensity from each state to change with the length of time spent in that state. The maximised likelihood changes from −1583 under the Markov model to −1481 under the semi-Markov model, while the estimated time spent in MDA over 10 years increases from 4.05 to 4.10. Given the estimates from this model, there is some evidence that the transition intensities, both to and from MDA, decrease with time spent in the current state. However a similar phase-type model with the partially-observed data would be challenging to define and identify from the data, and the principal results of interest appear to be robust to departures from the Markov assumption. The use of a multi-state model also allows a natural way to investigate the relationship between a composite outcome and explanatory variables. The model can be parameterised in terms of relative transition intensity functions for the multi-state model. However, it is also possible to make inference on relative measures of state occupancy which may provide a more useful presentation of the effects of explanatory variables. This approach to summarising findings from the use of a multi-state model may warrant consideration in other contexts.
As in Coates et al. (2010a) where the relationship between sustained MDA and the subsequent development of permanent damage in PsA was of interest, a composite outcome measure may also be of interest in terms of its longitudinal relationship to other outcomes. It is likely to be useful in this case to make use of a partially hidden multi-state model for the composite outcome, with its more comprehensive modelling, to understand better this relationship, particularly if prediction is not the only or primary focus of the investigation. In some cases, an approach to this might be to incorporate the partially hidden multi-state model framework into a larger multi-state model with state definitions also incorporating the additional outcomes to be related to the composite outcome. This has been done for the combination of simpler multi-state models in Tom and Farewell (2011). However, this will not always practical, or the most useful, approach, so further investigation of this problem is warranted.