All methods were carried out in accordance with relevant guidelines and regulations.

### Definitions and notation

Let states 1 and 2 represent the recovered and exacerbation states, respectively. Suppose *n* patients were followed over time to generate data on *n* independent alternating two-state processes. For patient *i* (\(i = 1, \ldots , n\)), let \(M_i\) be the number of observed exacerbations over their follow-up time \(T_i\). In our case study, the times \(T_i\) differ across patients because the nominal follow-up periods differed in the two trials under consideration, some patients dropped out during follow-up, and in both studies patients who were experiencing an exacerbation at the end of their nominal follow-up period were further followed to the termination of that event.

For the *j*th exacerbation of the *i*th patient, let \(U_{i,j}\) and \(V_{i,j}\) be, respectively, the time of exacerbation onset and termination. Thus, the exacerbation duration (within-exacerbation period) is \(W_{i,j}=V_{i,j}-U_{i,j}\) and the between-exacerbation period (the time from the termination of the previous exacerbation to the beginning of the *j*th exacerbation) is \(B_{i,j}=U_{i,j}-V_{i,j-1}\). In the recurrent event literature, these periods are often referred to as gap times. Since we do not know the termination time of the patient’s last exacerbation prior to study entry, we assume the randomization date coincides with the beginning of the initial exacerbation-free period; that is, we set \(V_{i,0} \equiv 0\). We will later discuss the implications of this simplifying assumption in the context of our case study. Similarly, we define \(S_{i,j}\) as the binary severity outcome of the *j*th exacerbation (\(S_{i,j}=1\) if the exacerbation is severe, \(S_{i,j}=0\) otherwise). We assume the censoring due to dropout occurs independently of the within- and between-exacerbation gap times and severity. In our case study, only the between-exacerbation gap times are censored as all individuals who were experiencing an exacerbation at the end of the nominal follow-up period were followed until the end of that exacerbation. Extension of our approach for situations with censored within-exacerbation gap times is immediate. Figure 1 provides a schematic illustration of the history of exacerbations in a given patient and the relationship between the three disease features of interest.

In addition, we let \(\varvec{X}_{i,j}\) (for \(j = 1, \ldots ,M_i\)) be the vector of covariates for the *j*th episode of patient *i*; the vector \(\varvec{X}_{i,j}\) includes both baseline covariates and episode-specific covariates that are measured at the onset of the *j*th exacerbation-free period. For simplicity of notation, in what follows we consider the same set of covariates \(\varvec{X}_{i,j}\) for modeling each of the three features, but optionally different sets of covariates can be used for each component.

### The model

We propose a parametric model to jointly characterize the distribution of \(B_{i,j}\) (between-exacerbation gap time, related to exacerbation rate), \(W_{i,j}\) (within exacerbation gap time or duration), and \(S_{i,j}\) (severity of exacerbation). Each of the three submodels incorporates a random effect to capture heterogeneity beyond that explained by included covariates, thus inducing autocorrelations among the repeated gap times of each type and the repeated severities within a patient. Joint modeling of the three random effects for each patient also allows for inter-correlations among the three model dimensions. We opted for parametric specification of all model components to allow straightforward estimation of random effects and prediction of outcomes. A further advantage of a parametric specification is that the likelihood function for the model can readily be programmed in available software.

#### The rate submodel

We use hazard-based models for the alternating two-state process [20]. The event hazard function gives the instantaneous rate of an event occurring at time *t*, conditional on the process history. For simplicity, we suppress the dependence on the process history and covariates in our notation. For the onset of an exacerbation event (the between-exacerbation gap times \(B_{i,j}\)), we use an accelerated failure time (AFT) model with an additive random effect (RE) [26, 27]. Specifically, conditional on a zero-mean individual-specific random effect \(Z_{B,i}\), the \(B_{i,j}\)’s of patient *i* (\(i=1, \ldots ,n\) and \(j=1, \ldots ,M_i\); for simplicity, the range of the indices will be suppressed hereafter) are assumed to be independent. The hazard function for \(B_{i,j}\) is modelled as

$$\begin{aligned} h_{i,j,B} \left( t | Z_{B,i} \right) =\theta _{i,j,B} h_B^0 \left( \theta _{i,j,B} \left( t-V_{i,j-1} \right) \right),\quad for\quad t > V_{i,j-1} \end{aligned}$$

where \(h_B^0 (.)\) is the baseline hazard function for the onset of exacerbation events and \(\theta _{i,j,B}=\exp \left( -\eta _{i,j,B} \right)\) where \(\eta _{i,j,B}=\varvec{X}_{i,j}^{T} \varvec{\beta }_B+ Z_{B,i}\), with \(\varvec{\beta }_B\) the covariate coefficient vector.

This parametrization of the AFT, which can equivalently be represented as

$$\begin{aligned} \log \left( B_{i,j} \right) = \eta _{i,j,B} + \epsilon _{i,j,B} \end{aligned}$$

where \(\exp \left( \epsilon _{i,j,B} \right)\) has hazard function \(h_B^0 (.)\), has the attractive feature that interpretations of the regression coefficients are the same either conditionally (on the random effect) or marginally [28]. We refer to this as an AFT-RE model.

Parametric time to event models require specification of the baseline hazard function. We suggest evaluating different functions and making a selection based on goodness-of-fit (e.g., as measured by the Akaike Information Criterion [AIC] [29]), as well as visual comparison of observed and predicted time to event curves. In contrast to hazard ratios estimated from proportional hazards models, the regression coefficients of AFT models can be interpreted intuitively and simply in terms of factors that contribute to the acceleration or deceleration of the time to events [30]. More precisely, with the parametrization of the AFT hazard function specified above, if \(\beta\) is the regression coefficient for a particular covariate, a one unit increase in that covariate is associated with multiplication of the expected time to event by the factor \(\exp (\beta )\). In what follows, we refer to these as AFT factors.

#### The duration submodel

We use another AFT-RE model for the termination of an exacerbation event (the within-exacerbation gap times \(W_{i,j}\)). Conditional on a zero-mean individual-specific random effect \(Z_{W,i}\), the \(W_{i,j}\)’s of patient *i* are assumed to be independent. The hazard function for \(W_{i,j}\) is modelled as

$$\begin{aligned} h_{i,j,W} \left( t | Z_{W,i} \right) = \theta _{i,j,W} h_W^0 \left( \theta _{i,j,W} \left( t-U_{i,j} \right) \right)\quad for\quad t > U_{i,j} \end{aligned}$$

where \(h_W^0 (.)\) is the baseline hazard function for the termination of exacerbation events and \(\theta _{i,j,W} = \exp \left( -\eta _{i,j,W} \right)\) where \(\eta _{i,j,W} = \varvec{X}_{i,j}^{T} \varvec{\beta }_W+ Z_{W,i}\), with \(\varvec{\beta }_W\) the covariate coefficient vector. This parametrization of the AFT can equivalently be represented as

$$\begin{aligned} \log \left( W_{i,j} \right) = \eta _{i,j,W} + \epsilon _{i,j,W} \end{aligned}$$

where \(\exp \left( \epsilon _{i,j,W} \right)\) has hazard function \(h_W^0 (.)\). Again, one can evaluate several parametric baseline hazard functions and make a selection based on goodness-of-fit.

#### The severity submodel

We specify a logistic regression model for the conditional probability of an exacerbation being severe once it occurs, \(p_{i,j} = P \left( S_{i,j} = 1 | Z_{S,i} \right)\), as

$$\begin{aligned} logit \left( p_{i,j} \right) = \varvec{X}_{i,j}^{T} \varvec{\beta }_S + Z_{S,i} \end{aligned}$$

where \(\varvec{\beta }_S\) is the covariate coefficient vector. Conditional on a zero-mean individual-specific random effect \(Z_{S,i}\), the \(S_{i,j}\)’s of patient *i* are assumed to be independent.

#### The zero-inflated component

In many circumstances, there is an excess number of patients who experience no exacerbations over the follow-up period (either because they are not susceptible to exacerbations or have a very low rate of events - a so-called “ non-susceptible” subgroup). This is similar to the situation motivating ‘cure models’ in survival analysis to model the presence of long-term survivors [31, 32]. To accommodate this aspect of the population, one can consider a mixture component which models \(\pi _i\), the probability patient *i* is non-susceptible, through a logistic regression model in terms of the baseline covariates \(\varvec{X}_{i,0}\) as

$$\begin{aligned} logit ( \pi _i ) = \varvec{X}_{i,0}^{T} \varvec{\beta }_{ZI} \end{aligned}$$

where \(\varvec{\beta }_{ZI}\) is the covariate coefficient vector. As there are no repeated events related to this mixture modeling, we do not include a RE in this model component.

#### Modeling the interdependencies

We assume for patient *i*, given the random effects, \(B_{i,j}\) and (\(W_{i,j},S_{i,j}\)) are independent, and given the random effects and \(B_{i,j}\), \(W_{i,j}\) and \(S_{i,j}\) are independent (for \(j=1, \ldots , M_i\)). These assumptions allow for two types of correlation between the outcomes based on the random effects. The random effects in each submodel account for autocorrelation in that outcome across repeated events on the same patient. In addition, to accommodate potential dependencies among the three outcomes (e.g., if individuals with higher exacerbation rates tend to experience more severe exacerbations), we allow for correlation between the random effects across submodels. Specifically, we assume \(\varvec{Z}_i = (Z_{B,i}, Z_{W,i}, Z_{S,i} ) ^ T\), \(i=1, \ldots , n\), are independent and identically distributed with mean zero. For our motivating case study, we will take this distribution to be multivariate normal with covariance matrix \(\Sigma _{\varvec{Z}}\). The distribution of the random effects is thus governed by 3 variance and 3 correlation parameters; these correlation parameters indirectly describe the relationships across patients between the rate, duration, and severity of the exacerbations.

### Implementation

The likelihood function is provided in Section 2 of the Additional file 1. We use PROC NLMIXED in SAS (with the built-in Newton–Raphson ridge optimization algorithm) to obtain the maximizer of the full likelihood, together with its estimated variance-covariance matrix (using SAS 9.4 PROC NLMIXED [SAS Institute, Cary NC]). To obtain the AFT factors/odds ratios (ORs) and their confidence intervals (CIs), we exponentiate the estimated coefficients and their CIs from the linear predictors. An annotated SAS macro that implements the model in a generic fashion, along with a manual and a simulated dataset, is available at our website (http://resp.core.ubc.ca/software/RDSmodel). The macro flexibly accommodates various aspects of the model such as the inclusion or exclusion of the duration and severity submodels, or which submodels should include a random effect.