A multistate model for events defined by prolonged observation.

Time-to-event and similar analyses can be problematic if the event of interest is operationally defined by some condition being true for a prolonged period of time. A particular example of this, remission in psoriatic arthritis, is considered in detail for illustration. A 3-state model is proposed for characterizing the transition rates into and out of remission. Remission is linked to an initial and subsequent state for the purpose of introducing the condition that remission must be of some duration to be clinically meaningful. The model is compared with alternative approaches that have been used in such situations. These involve 2-state models where the duration of remission is allowed for through different definitions for the time of entry into remission. Both definitions are linked to prolonged observation of a particular clinical state.

INTRODUCTION

In some circumstances, a time-to-event analysis is desired for an event that is difficult to define observationally. In this paper, we will focus on a particular example motivated by a published study of disease remission for patients with psoriatic arthritis (PsA) (Gladman and others, 2001). In this study, patients were followed in a hospital clinic and were defined to be in remission if they were observed to have no evidence of disease activity for at least 3 consecutive clinic visits. Our aim was to illustrate an approach to the analysis of such data and the potential bias associated with the adoption of simple but possibly inappropriate definitions of the time to an event.

THE PSA DATA

We will make use of a clinical database from the University of Toronto Psoriatic Arthritis Clinic. The database contains information on 790 patients entering the clinic in the years 1973–2006. PsA is an inflammatory arthritis associated with psoriasis. A basic measure of disease activity is a clinically obtained count of the number of joints that are swollen and/or painful, that is the active joint count. These counts may vary considerably across clinic visits and zero counts are not uncommon. Thus, from a clinical perspective, it is only a long period of inactive disease that can be regarded as time in remission. There is interest in factors associated with the chance of remission as information on these may contribute to treatment decisions.

Rather than define an observed time to remission, we choose to regard remission as a conceptual disease state with those patients not in remission assumed to be in an active disease state. Observationally however, we link remission to the active joint count in the following manner. If a patient has 1 or more active joints at some time point, then their disease state is active. Thus, any patient in remission must have no active joints and this can be regarded as a structural zero. If a patient has had 3 consecutive clinic visits when no active joints were observed, then the patient is regarded as being in remission at the last clinic visit during this period. Note that this consistent observation is not taken to define remission but is taken to “confirm” that the patient has entered remission at some point prior to the third visit. Given the intended pattern of clinic visits, 3 consecutive clinic visits will usually take place over a period of at least a year and will almost always be longer than 6 months since usually there must be at least 3 months between visits in order for the second to be regarded as a new visit requiring a protocol to be completed. It would be possible to require that zero counts be observed at all clinic visits over some defined time interval, or greater interval, be taken as confirmatory but it is slightly simpler to link confirmation to a fixed number of visits which, for these data, correspond to a minimal time period at least.

A patient with no active joints at a clinic visit that is not preceded by at least 2 other such visits may be in remission, but it cannot be confirmed. They may also have active disease with the observed zero joint count being a sampling zero. It is well recognized clinically that sporadic zero joint counts can occur for patients with active disease. We will incorporate this uncertainty in the approach introduced in Section 3.

In practice, 2 other approaches can be used to define remission in the PsA data. The first treats a patient as in remission at the end of a time period with 3 consecutive clinic visits with no active joints, while in the second approach a patient was considered in remission at the beginning of such a period. We will compare all 3 approaches in the illustration in Section 4.

A STATISTICAL MODEL

A simple 2-state reversible multistate model in continuous time, with one state corresponding to active disease (i.e. with observed active joints) and the other to remission (i.e. without active joints), would provide a structure for the analysis of time to remission. However, without careful modeling of the backward transition from remission to active disease, such a model will fail to capture the essential characteristic of remission, that is the time in this state cannot be of a short duration if it is to have clinical meaning. An assumption that transition out of remission is impossible prior to some fixed time in the state, say 6 months, could, in principle, be made. However, modeling such a transition rate would require special software and the assumption of any single minimal time period for all patients would be somewhat arbitrary.

To overcome this limitation of the usual 2-state model, we propose a 3-state model as illustrated in Figure 1. Suppose that S(t) represents the state for a patient at time t, where t denotes the years since diagnosis of PsA. In this model, state 1 (S(t) = S1) corresponds to active disease, state 2 (S(t) = S2) corresponds to the “early” stage of remission and state 3 (S(t) = S3) corresponds to an “established” stage of remission. It should be noted that a patient in either S2 or S3 is regarded as in remission. These are conceptual, not observed, states that are introduced in order to give all remissions some duration. This is achieved by the further assumption that patients in S2 can only move to S3. Transitions back to active disease, S1, can only occur from S3. This model can be further specified by 3 transition rates or hazard functions. Associated with the S1→S2 transition, the first is the instantaneous hazard of progression to S2, conditional on occupying S1 at t:where λ12,0(t) is a baseline hazard, x represents a vector of explanatory variables and β12 is the corresponding unknown vector of regression coefficients. In the example under consideration here, this transition into S2 would be the transition of greatest interest as it would provide a model for the distribution of time to remission from an active disease state.

Fig. 1.

Three-state model for remission in PsA1

A hazard for the S3→S1 transition, λ31(t,x) can be defined analogously and corresponds to a model for transitions from remission back to active disease. The final transition, S2→S3, will have a hazard of the same form but this conceptual transition is unlikely to be of clinical interest. In principle, its estimation provides information on the duration of remission, along with λ31(t,x), but S2 is used here primarily for convenience. If the modeling of remission duration was of direct interest, then alternative approaches might be of more value (Temkin, 1978).

The conceptual nature of this model requires that uncertainty regarding a patient's state at some time point be made explicit and this can be done through the introduction of misclassification probabilities. Suppose that S(t) = r (r = S1,S2,S3) represents the true underlying state of a patient at time t, and O(t) = s (s = O1,O2,O3) represents the corresponding observed state at time t based on active joint counts. Specifically, let O1 denote the observation of a nonzero active joint count at a visit, O2 denote the observation of a zero active joint count at a visit not preceded by at least 2 other such visits, and O3 denote the observation of the third or subsequent zero active joint counts in a sequence of such visits. Then, based on the assumptions about active joint counts outlined in Section 2 and the assumption that the misclassification probabilities are independent of time t, we can specify Pr{O(t) = s∣S(t) = r} as follows: (1)

(1) Pr(O1∣S1) = 1 − Pr(O2∣S1), Pr(O2∣S1) = logit − 1(γz), Pr(O3∣S1) = 0, where z is a vector of explanatory variables that might be informative about the true state of the patient at time t, and γ is the corresponding unknown vector of regression coefficients; (2)

(2) Pr(O1∣S2) = 0, Pr(O2∣S2) = 1, Pr(O3∣S2) = 0; (3)

(3) Pr(O1∣S3) = 0, Pr(O2∣S3) = 0, Pr(O3∣S3) = 1.

Essentially, in this model, we only allow S1 to be misclassified. In other words, at some time point t, a patient with a clinic visit and an associated zero active joint count not preceded by at least 2 other such visits can either be in active disease (S1) or be in early stage of remission (S2). Furthermore, we use a logistic model to investigate the influence of explanatory variables z on Pr(O2∣S1). Specifically, for the illustration in Section 4, we will include a binary explanatory variable, Z, that is coded 1 for the first visit with a zero joint count not preceded by any other such visits and 0 for visits with a zero active joint count only preceded by one such visit. The variable is defined only for visits corresponding to observed state O2. Other information such as the gaps between visits could be considered but, for simplicity, this is not investigated here.

Multistate models with misclassification have been discussed in Jackson and others (2003). The estimation of misclassification probabilities derives from the joint estimation of the model for misclassification and the transtion rates for the multistate model for the true states. Misclassification is not objectively determined for any patient but patients who have 1 or 2 zero counts, corresponding to O2 observations, are potentially misclassified and will contribute terms involving misclassification to the likelihood. As outlined in Jackson and others (2003), the likelihood contribution from an observed state O is ΣPr(O∣S)Pr(S),j = 1,2,3, with the likelihood contribution from an individual being the product of such terms for each observation time. Heuristically, information on misclassification thus derives from arguments such as: if a patient has a single zero observation and is in remission, then the patient must go through states S2 and S3 before returning to S1 by the next clinic visit. If the time period between visits is relatively short, then the zero observation is unlikely to derive from the remission state.

These models can be fit using the msm package in R (http://www.r-project.org) under the assumption that transition rates are constant (i.e. time independent) or piecewise constant over time. For convenience therefore, we will illustrate the use of our remission model with this assumption. Because we only observe the (O1,O2,O3) classification at clinic visits, the maximum likelihood estimation is based on the panel data arising from this intermittent observation pattern. However, interval censoring of transition times is incorporated as an option in the msm package.

ILLUSTRATION

Models

We fit the 3-state model with misclassification (denoted Model A) to the Toronto PsA data introduced in Section 2. The time scale t in the model is the years since diagnosis of PsA, that is the disease duration, which is typically > 0 at clinic entry because patients may be referred or come to the clinic some time after their diagnosis. This left truncation is easily handled for constant or piecewise constant transition rates and can, in principle, be handled for more general transition rates in multistate models (e.g. Siannis and others 2007). For S1→S2 and S3→S1 transition rates, the explanatory variables included are sex (coded as 1 for males, 0 for females) and age at PsA diagnosis (standardized by the sample mean 35 years and sample standard deviation 12 years). No explanatory variables are included for the transition S2→S3. As indicated earlier, for the logistic model of misclassification probability, a binary explanatory variable Z coding the first visit with zero active joint count not preceded by other such visits is used. It is assumed that at t = 0, all patients were in S1, the active disease state; and all transitions are interval-censored due to the intermittent observation pattern of clinic visits. For simplicity, we assume that the interval censoring and any right censoring due to mortality are noninformative in this illustrative example.

Preliminary analyses indicated that the assumption of constant baseline transition rates over time might not be appropriate as the transitions into remission were more likely to occur in later follow-up periods, given the explanatory variables. To account for this phenomenon, we allow the baseline transition rates to be piecewise constant such that t = 15 (the 60th percentile of observed disease durations) is a change point for these rates. Models fit with change points at 5, 10, 15, and 20 years indicated that the significant change was at 15 years and hence a simple model with this single change point is used for illustration. Methods for the estimation of change points in exponential models have been investigated (Matthews and others, 1985) but the piecewise constant assumption in multistate models is usually adopted as a simple representation of a nonconstant hazard and not because there is particular inferential interest in a change point.

Two other multistate models with only 2 states, corresponding to S1 and S3 in the 3-state model, and no misclassification (denoted Models B and C) are also fit to the PsA data. In Model B, it is assumed that patients were in remission only at the third of 3 or more consecutive clinic visits with zero active joint counts (either from enrollment or following a visit with nonzero active joint count). That is, a transition to remission had occurred between the second and third clinic visits of a time period with 3 consecutive visits with zero active joint counts. In Model C, it is assumed that remission had occurred by the first visit of such a period and after the previous visit at which a nonzero count was observed. These represent modeling strategies that might be adopted for these remission data in practice. Again, for the rates associated with the S1→S3 and S3→S1 transitions, sex and age at PsA diagnosis are included as explanatory variables; and we allow baseline transition rates to be piecewise constant with a change point at t = 15.

To illustrate the different assumptions made by the models, Figure 2 displays the time periods where remission may be assumed to have started for the 3 models, for the situation when a positive active joint count is followed by 3 consecutive zero counts.

Fig. 2.

Time periods where remission can start in Models A, B, and C in the scenario when a visit with a positive active joint count was followed by 3 consecutive visits with zero joint counts. 2

For Models B and C, there is a choice involved in the handling of individuals who enter the clinic with a zero joint count, and those whose last 2 visits have had zero joint counts preceded by a nonzero count. Excluding initial visits with a zero joint count leads to little change in model estimation and these have therefore been included. One or 2 visits at the end of follow-up would not be taken to represent remission in application of these approaches and therefore these have also been included. Inclusion of all these visits ensures that the same data are being used for all model fits.

Results

Observed transitions.

The frequencies of observed transitions between the states in Model A, B and C are given in Table 1. Furthermore, in Model A, there were 443 observations of 2 consecutive visits with zero active joint counts from enrollment or following a visit of nonzero active joint count; 255 of them were further followed by a third visit with zero active joint count, 134 of them were followed by a visit with nonzero active joint count, and 54 of them were observed at the last 2 clinic visits up to the current study termination date.

Table 1.

Frequencies of observed transitions between states in Models A, B and C

Model A			To state
			1	2	3
	From state	1	4642	890	0
		2	514	443	255
		3	168	0	447
Model B			To state
			1	3
	From state	1	6489	255
		3	168	447
Model C			To state
			1	3
	From state	1	5994	240
		3	168	957

Estimated transition rates.

The estimated baseline transition rates and effects of explanatory variables from the 3 models are given in Table 2. Because the estimated effects of explanatory variables vary across models that additionally have different underlying assumptions, comparison of the estimated transition rates is not particularly informative. Nevertheless, it can be noted that, for all 3 models, the estimated transition rates both to remission and back to active disease were higher in the later time period (t∈[15,∞)). With respect to the estimated regression coefficients associated with sex, these are similar for the 3 model fits: given the same time period and age at PsA diagnosis, male patients were more likely to have transitions to remission and less likely to have transitions back to active disease. A similar sex effect was also reported in Gladman and others (2001). In contrast, an effect of age at PsA diagnosis was not evident in Models B and C, while in Model A, older patients at PsA diagnosis tended to have more transitions both to remission and back to active disease. There is some evidence that older patients have more observations of 1 or 2 visits with a zero count followed by a visit with a positive count, and that long gap times occur less frequently for these patients. These features may be linked to the differential results concerning age but this is difficult to ascertain for certain and it seems best to regard the differential results as simply raising a question about the robustness of the lack of an age at diagnosis effect observed in Models B and C and in Gladman and others (2001).

Table 2.

Parameter estimates, followed by 95% confidence intervals in parentheses, from Models A, B, and C

Model A
	Baseline transition rates		Time period: [15, ∞]	Sex	Age at PsA diagnosis
	λ12, 0: 0.042 (0.033,0.053)		0.283 (0.010, 0.555)	0.614 (0.344,0.884)	0.177 (0.048,0.306)
	λ31, 0: 0.328 (0.225,0.478)		1.375 (0.950, 1.800)	-0.198 (-0.524,0.128)	0.458 (0.295,0.622)
	λ23, 0: 1.118 (0.954,1.311)		-0.921 (-1.243,-0.598)
	Misclassification
	Pr(O2 | S1)	0.031 (0.065,0.087)
	γ	8.069 (7.356,8.783)
	Total length of stay for t ∈ [0, 40]
			State 1	State 2	State 3
	Male, age at PsA diagnosis = 35		30.593	5.140	4.267
	Female, age at PsA diagnosis = 35		34.647	3.157	2.196
Model B
	Baseline transition rates		Time period: [15, ∞]	Sex	Age at PsA diagnosis
	λ13, 0: 0.041 (0.028,0.059)		0.603 (0.324,0.881)	0.747 (0.449,1.045)	0.058 (-0.085,0.201)
	λ31, 0: 0.696 (0.473,1.024)		0.071 (-0.252,0.394)	-0.318 (-0.638,0.003)	-0.074 (-0.238,0.089)
	Total length of stay for t ∈ [0, 40]
			State 1		State 3
	Male, age at PsA diagnosis = 35		32.467		7.533
	Female, age at PsA diagnosis = 35		36.981		3.019
Model C
	period: [15, ∞]		Sex	Age at PsA diagnosis
	λ13, 0: 0.034 (0.024,0.049)		0.614 (0.340,0.888)	0.830 (0.538,1.122)	0.102 (-0.038,0.242)
	λ31, 0: 0.300 (0.203,0.442)		0.320 (0.001,0.638)	-0.219 (-0.540,0.102)	0.075 (-0.086,0.236)
	Total length of stay for t ∈ [0, 40]
			State 1		State 3
	Male, age at PsA diagnosis = 35		29.482		10.518
	Female, age at PsA diagnosis = 35		35.451		4.549

Total length of stay in active disease and in remission.

A more revealing comparison of the fitted models is when they are used to estimate the total length of stay in the various states over a fixed period (Table 2). The calculation of these quantities has been discussed in Jackson (2007). Following van den Hout and Matthews (2008), we incorporate the change in baseline transition rates over the time period (t∈[0,15) and t∈[15,∞)) for these calculations.

If we choose a period from t = 0 to t = 40 (the 98th percentile of the observed disease durations) and assume no mortality, then, for male patients who were 35 years old at PsA diagnosis, the estimated time with active disease is 30.6, 32.5 and 29.5 years in Models A, B, and C, respectively. Correspondingly, the time in remission (time spent in either S2 or S3 in Model A) for the same patients follows the reverse ordering. The results for female patients who were 35 years old at PsA diagnosis are different: the estimated time with active disease is 34.6, 37.0, and 35.5 years in Models A, B, and C, respectively.

It can be seen that as the most conservative approach of defining remission for these data, Model B always gives the least time spent in remission, irrespective of explanatory variables. Model A utilizes more information from the visits with zero active joint counts but not preceded by at least 2 other such visits, and allows the possibility that at these visits the patients were actually in early stage of remission. Therefore, the estimated time spent in remission in Model A is expected to be longer than in Model B. In Model C, we allow the patients to be in remission exactly 2 clinic visits earlier than in Model B, which also makes the estimated time spent in remission longer. The ranking of the corresponding estimates from Model A and Model C will depend on specific scenarios. In Model A, we use the information of the visits with zero active joint counts but not preceded by at least 2 other such visits, regardless of the observed states at the following visits. In Model C, given that there are at least 3 consecutive visits with zero active joint counts, we allow the patients to be in remission as early as at the first visit of such a sequence of visits (without consideration for misclassification). However, we do not consider the possibility of remission at those 1 or 2 visits with zero active joint counts not preceded by at least 2 other such visits and immediately being followed by a visit with nonzero active joint count or reaching the study cutoff date. There exists a trade-off between these factors, which determines the actual estimates from Model A and Model C for the length of stay in the various states.

Misclassification probabilities.

It is also of interest to examine the estimated misclassification probabilities in Model A. Compared with Pr{O(t) = O2∣S(t) = S1,Z} parameterized in Model A, a more clinically relevant measure would be Pr{S(t) = S1∣O(t) = O2,x,Z}, that is the probability that patients were actually in active disease when observed with clinic visits with zero active joint counts not preceded by at least 2 other such visits, given the explanatory variables. This can be obtained through Bayes' rule:

Since Pr{O(t) = O2∣S(t) = S3,Z} = 0, Pr{S(t) = S1∣x}, and Pr{S(t) = S2∣x} are needed for the above calculation and they can be approximated by the proportions of estimated total length of stay in S1 and S2 for the time periods t∈[0,15) and t∈[15,40], given the explanatory variables. For females diagnosed at age 35, these proportions were (0.88, 0.03, 0.09) in the first time period and (0.85, 0.11, 0.04) for the second. For males, the corresponding proportions were (0.79, 0.05, 0.16) and (0.75, 0.18, 0.07). The resultant missclassification probabilities are shown in Table 3.

Table 3.

Approximate misclassification probabilities Pr{S(t) = S1|O(t) = O2} by explanatory variables in Model A

Misclassification probability Pr{S(t) = S1|O(t) = O2}			First zero	Second zero
Age at PsA diagnosis = 35	Male	Time period: [0, 15)	0.938	0.322
		Time period: [15, 40]	0.810	0.118
	Female	Time period: [0, 15)	0.965	0.468
		Time period: [15, 40]	0.887	0.199

It can be seen that at the second of the visits with zero active joint counts not preceded by at least 2 other such visits, the misclassification probability Pr{S(t) = S1∣O(t) = O2} is much smaller than at the first of such visits. In addition, misclassification was less likely to occur for t∈[15,40], which is consistent with our findings of higher transition rates to remission in the later follow-up period. Male patients were also less likely to be misclassified, which is expected since male patients also tended to have higher transition rates to remission than females. Note that in the model for Pr{O(t) = O2∣S(t) = S1,Z}, the effect of the binary explanatory variable Z is large compared with the one evident for Pr{S(t) = S1∣O(t) = O2,x,Z}. This is due to the large Pr{S(t) = S1∣x} and small Pr{S(t) = S2∣x} estimated from the PsA data. That is, the patients primarily had active disease during the follow-up period.

DISCUSSION

There are a variety of contexts in which interest lies in the distribution of times to events which are manifested by a defined period of time in a particular state. In addition to the example of remission in PsA, a second example, that also motivated the work reported here, concerned an analysis of times to being lost to follow-up for patients in a rheumatoid arthritis registry (Krishnan and others, 2004). An accompanying editorial (Farewell, 2004) pointed out the potential difficulties with the usual definitions adopted for the actual time at which lost to follow-up occurred. The difficulties arise, because, either explicitly or implicitly, the state of being lost to follow-up will usually be defined by nonacquisition of new data for a period of time.

We have considered a 3-state model (Model A) for the analysis of the time to events which are manifested by a defined period of time in a particular state. The initial time observed in the state (Model C) or the end of the defined period (Model B) could be regarded as the observed time to the event (or the first time a patient is observed to have had the event with panel data). Both these strategies have drawbacks. The use of the 3-state model avoids these drawbacks but does introduce the need for a somewhat “hypothetical” transition after entering the state of interest in order to allow more plausible times in the state. However, if misclassification probabilities are included in the model, then the model avoids the need, required under other strategies, to assign observed states that are known to be problematic.

Our consideration of other strategies has been based on a panel data assumption. The assumption that clinic visits correspond to exact transition times when fitting Models B and C would not change the qualitative conclusions. An additional alternative would be to fit Model B under the assumption that a patient enters remission between their last observed visit with a nonzero joint count and the third visit with an observed zero joint count. When used with the PsA data, this gives observed times in states comparable to our 3-state model and explanatory variable effects similar to the Model B analysis. This approach, as in Model B, does not allow any patient to be potentially regarded as in remission unless 3 consecutive visits with zero joint counts are observed.

Minimally, the use of the 3-state model could be seen as a check on the robustness of findings from more simplistic approaches that might be adequate for some purposes. Routine use of the methodology will require computational tools to fit the models with more complicated or possibly semiparametric hazard models.

FUNDING

Medical Research Council, UK (U.1052.00.009).