Depending on the research question, different approaches can be used for the statistical analysis of antidepressant medication data. If the focus of interest are rates until first medication use (medication initiation), a single event survival analysis suffices. However, for other research questions, such as the probability of being in a medication cycle since the start of the follow-up or since entering a medication discontinuation state, specific multi-state structures should be used. Figure 1 provides an overview of the multi-state structures used in the current study, from the simplest structural approach (Fig. 1A Single event survival analysis) to the most complex model (Fig. 1G: MSM with recurrent medication cycles/medication discontinuation periods with restrictions applied). The graphs presented in Fig. 1 were produced via the interactive web-tool MSMplus (23). The traits, results, interpretation, advantages and drawbacks of each structure are presented in section 4.4.

A) Single-event survival analysis, B) Competing risks, C) 3 –state Illness-Death model, D) Multi-state model with medication discontinuation state (4-state Unidirectional model), E) Bidirectional structure between medication cycles and discontinuation periods (4-state Bidirectional model), F) MSM structure with recurrent medication cycles and discontinuation periods, G) Recurrent MSM structure with restrictions - “Emulated bidirectional” structure.

The list above is not an exhaustive list of multi-state structures that can be used and research questions that can be addressed. The estimated measures can be compared between the exposure groups of interest.

## 4.4 Structures application, results and interpretation

**Single event survival analysis**

When the interest lies in time until a single event of interest (e.g time to first antidepressant medication), a simple hazard model suffices. In this case, the individual is at risk of transitioning from an initial, medication-free state to the first medication cycle (medication initiation) (Fig. 1A). In this type of analysis, the rate of that transition is the estimate of interest, as well as the relevant risk differences and risk ratios between different groups. Figure 2a shows the estimated survival probabilities (not having medication up until time \(t\)) for the group of matched healthy-controls and the BC cases while Fig. 2c shows the probability of medication initiation up to time \(t\)which is \(1-S\left(t\right)\). The cases have a higher probability for antidepressant medication initiation, with a 20 % probablity for at least one prescription of antidepressants compared to 10 % for th matched healthy controls within 6 years since start of follow-up. These are the probabilities of ever been prescribed medication in the hypothetical situation that the individual cannot die due to any causes and thus tend to overestimate the true probabilities of ever been prescribed medication. Figure 2b show that BC cases have almost 3 times higher hazard rate for first medication use at 1 year since start of follow-up compared to the healthy-matched individuals.

Medication-free survival probability (a), Hazard ratio for antidepressant medication initiation (b), probability of antidepressant medication initiation (c) for controls and cases.

The single-event survival analysis approach is straightforward, easy to implement, and fine-tailored for comparison of rates. However, it fails to account for the fact that individuals diagnosed with BC have higher risk of dying compared to healthy individuals. Due to this fact, the derived probabilities potentially overestimate the outcome of interest, that is, the probability of ever antidepressant medication initiation up until time \(t\). The approach that follows considers the competing risk of death when estimating probabilities of antidepressants initiation.

Competing risks

In a competing events structure, the individual found in the initial medication-free state can experience either the medication initiation state or the competing event of death (Fig. 1B). This approach takes into account the competing risk of death in the estimation of the probability of experiencing the first antidepressants medication. We can interpret this probability as “the probability of antidepressant medication initiation up to time \(t\) after the start of the follow-up”. In Fig. 3a, the probability of ever been prescribed medication (medication initiation) is derived for the cases and the controls, having accounted for the fact that the rate of death is much higher for the women diagnosed with BC (red versus blue dash-dot lines of Fig. 3c). It can be observed that in this setting, the estimates from the single-event survival analysis are very close to the ones derived by the competing risks approach (Fig. 2c versus Fig. 3a), as death due to BC does not seem to be a strong competing event for time until antidepressant medication initiation. Even though the competing risks approach considers the competing risk of death, it does not permit the individual to leave the state of experiencing the first medication use. Thus, more complex multi-state structures are in need if the interest lies in what happens after the first medication use.

3-state Illness- Death multi-state model

If a transition between the medication initiation and the absorbing state of death is allowed, the individual can experience death after medication initiation (Fig. 1C). This structure allows the study of the probability of having experienced the first medication use while also allowing individuals who have entered this state to move to the death state. Thus, the derived probability in Fig. 3b can be interpreted as the probability of ever been prescribed medication and still be alive up to time \(t\) after the start of the follow-up. For the controls, the transition rate from medication initiation towards death is low (Fig. 3d) so the probability estimates of Fig. 3b are very close to those of Fig. 3a (Competing risks structure). However, for women diagnosed with BC, the transition rate towards death both before and after the first medication use is higher (Fig. 3d). For this reason, the estimated probabilities of experiencing the first medication use and still being alive for BC cases are lower than the probability of experiencing the first medication use derived from the competing risks approach.

The hazard rates towards death for the controls (with or without medication) are approximately 0.001

The medication initiation state consists of people who either stayed in a medication cycle for the rest of their follow-up or, more likely, those who moved on to discontinuation periods or subsequent medication cycles. However, the 3-state Illness-Death multi-state model is not able to discern this issue and thus is not appropriate for studying the probability of being in a medication cycle and the length of stay in that cycle. A multi-state model that includes a discontinuation period state would allow for the study of being in the first medication cycle.

Multi-state model with a medication discontinuation period state

By adding a discontinuation period-state (Fig. 1D), the individual is allowed to either go to the absorbing state of death or to a medication discontinuation period after entering the 1st medication cycle state. This addition allows us to shift our focus from the probability of ever been prescribed medication to “the probability of being in the 1st medication cycle” at time \(t\) since the start of follow-up or time since entering the 1st medication cycle. Subfigure 4a depicts the probability of being in the 1st medication cycle as a function of time\(t\) since the start of follow-up while subfigure 4b shows the probability of being in the 1st medication cycle and the probability of heading towards the 1st medication discontinuation period as a function of time since entering the cycle. The overall probability of being in the 1st medication cycle (Fig. 4a) is higher for the cases (red) compared with the healthy-matched individuals (blue line). The probability of staying in the 1st medication cycles state is higher for the cases for the first 1.5 years since entering the medication cycle.

Probability of being in first medication cycle since the start of the follow-up (a), Probability of staying (solid) in the first medication cycle and entering the 1st discontinuation period (dash) for time t since entering the cycle (b).

Under this structure, the estimated probability for the 1st medication discontinuation period can be interpreted as “the probability of ever having exited the 1st medication cycle and still be alive up to time \(t\)”, consisting of people who either stayed in a medication discontinuation state for the rest of their follow-up or experienced additional medication cycles later on during follow-up. Thus, if probability of being in the 1st medication discontinuation period and subsequent medication cycles and medication discontinuation periods is of interest, we need to allow an individual to be able to leave the 1st discontinuation period and enter in a new medication cycle.

Bidirectional model

A bidirectional multi-state structure can be built by allowing an individual to re-enter a medication cycle after they enter a medication discontinuation period (Fig. 1E). Therefore, a back transition is inserted from the discontinuation period state to the medication cycle state, allowing the individual to move back and forth multiple times -without a limit- among these states. The existence of the fourth state (discontinuation period) and the double transition arrow between this state and the medication cycle state serves a certain purpose. Allowing for a medication-free state that is different from the initial “Start of follow up” state and putting a back transition towards a medication cycle allows the transition rate for a next medication cycle to differ from the transition rate from the starting state to the first medication cycle. Figure 5 depicts the probability of being in a medication cycle since the start of the follow-up (5a), and the probability of being in a medication cycle or discontinuation period, given that an individual starts in a medication cycle(Fig. 5b) or a discontinuation period state (Fig. 5c). In Fig. 5a, the BC cases present a higher overall probability of being in a medication cycle compared to the controls over the follow-up time. In Fig. 5b it can be observed that BC cases have higher probabilities of being in a medication cycle and lower probabilities of being in a medication discontinuation period compared to controls as a function of time since entering a medication cycle. From Fig. 5c, starting in a medication discontinuation period, the BC cases are less likely to stay in a discontinuation period and more likely to be in a medication cycle compared to controls as a function of time since entering a discontinuation period.

Probability of being in a medication cycle at time t since the start of the follow-up (a), Probability of being in a medication cycle (solid) or being in a medication discontinuation period (dash) as a function of time t since entering a medication cycle (b), Probability of being in a medication cycle (dash) or being in a medication discontinuation period (solid) as a function of time t since entering a medication discontinuation period (c).

Despite its simple application and interpretation, this structure imposes the same transition rates between a medication cycle and a discontinuation period (and vice-versa) irrespectively of the number of previous medication cycles (or discontinuation periods). While a time-varying covariate of previous medication cycles could be added in the transition rate models, the estimation of measures other than the transition rates and ratios is not possible under the Markov or the Semi-Markov assumption. This means that an individual has the same transition rate from a medication cycle to a discontinuation period (and vice-versa) no matter how many previous cycles they may have experienced, an assumption that is not very realistic in our setting. A recurrent multi-state structure allows for the separate modelling of a transition rate between each subsequent couple of medication cycle and discontinuation period (and vice-versa) thus accounting for past transitions via its own structure and is described below.

Recurrent multi-state model

A flexible approach when dealing with recurring states is to fit a recurrent multi-state structure that consists of repeated, ordered events/couples of medication cycles- discontinuation periods (Fig. 1F). This structure can be more flexible than the bidirectional model structure, as it allows the separate modelling of each transition rate without imposing same transition rates every time an individual moves from a medication cycle to a discontinuation period (and vice-versa). Due to recurrent states gradually becoming sparsely populated, there is a need for setting a threshold to the number of medication cycles and discontinuation periods to be modelled. In the case of this study, a maximum of six medication cycles are modelled, with subsequent medication cycles and discontinuation periods being ignored. Individuals that enter the 6th medication cycle are characterized as chronic antidepressant medication users and stay in this state until the end of their follow-up or move towards the death state.

By applying the current multi-state structure, there are estimated probabilities with useful clinical interpretation such as “What is probability of being in a medication period or discontinuation period since the start of follow-up (total probability on a population level) or given entering the \({1}^{st}, {2}^{nd},\dots , {k}^{th}\)medication cycle or discontinuation period?” Fig. 6a depicts the total probability of being in a medication cycle/the prevalence of medication use in the population since the start of follow-up on a population level. While this probability in absolute terms is low, cases appear to have more than three times the probability of being in a medication cycle since their diagnosis compared with controls for the first 2 years of their follow-up. Figure 6b depicts the probabilities for cases and controls of staying in their \({1}^{st}\), \({2}^{nd}\) and \({3}^{rd}\) medication cycle after entering each one of them. The probability of staying in a medication cycle increases with entering each new medication cycle for both cases and controls, with cases having higher probabilities compared to controls in all the three cycles presented here. Figures 6c and 6d present the total probability of a participant being in a medication cycle (both the current one and all the subsequent ones) as a function of time since entering the \({1}^{st}\), \({2}^{nd}\) and \({3}^{rd}\) medication cycle (Fig. 6c) or the \({1}^{st}\), \({2}^{nd}\) and \({3}^{rd}\) discontinuation period (Fig. 6d). The total probability of being in a medication period seems to increase as the number of past medication cycles increases, both as a function of time since entering a medication cycle or time since entering a discontinuation period. An interesting finding is that, given that a participant experiences her 3rd medication cycle or discontinuation period, a control has higher total probability of being in a medication cycle compared to a case.

Total probability on a population level of being in a medication cycle since the start of follow-up (a), Probability of staying in current medication cycle since entering it (b), Total probability of being in a medication cycle since entering the 1st, 2nd, 3rd medication cycle (c) or the 1st, 2nd, 3rd medication discontinuation period (d)

Using this approach allows the derivation of clinically significant measures (probabilities and length of stay) taking into account past medication cycles via the multi-state structure itself, a trait that the aforementioned bidirectional model lacked. However, using recurrent multi-state structures comes with its own drawbacks. A small to moderate increase in the number of states can greatly increase the structural complexity. Apart from the starting state (start of follow-up) and the absorbing state (death), there will be K medication cycle states (\({1}^{st}, {2}^{nd},{3}^{rd},\dots , {k}^{th}\)) and \(K-1\) post-medication period states (\({1}^{st}, {2}^{nd},{3}^{rd},\dots , {k}^{th}-1\)), leading to transitions having sparse data as state order progresses. This data sparsity can make transition-specific estimations troublesome, leading to convergence issues and low precision. Applying restrictions to the transition rates via shared parameter estimation for certain transitions, can address, at least partially, these issues. Keeping the modelling of the covariate effects and the baseline transition rates simple can help dealing with these issues. Applying restrictions to the transition rates via shared parameter estimation for certain transitions, can also address, at least partially, these issues.

Emulating a bidirectional structure by imposing restrictions

Under reasonable assumptions, sharing information across transitions of a multi-state structure in order to address the data sparsity issue described above is possible by imposing constraints in the parameter estimation. For the recurrent multi-state structure described (Fig. 1F), it can be assumed that transition rates towards death are not influenced by the number of past medication cycle states. This translates in transition rates from all medication cycle states to death being restricted to be the same and transition rates from all discontinuation period states to death to also be the same. Additionally, it can be assumed that each new transition between medication cycle and discontinuation period (and vice-versa) has a common underlying relationship with the time to experiencing the transition event (same shape) but on different scale (proportional). This translates to imposing commonly shaped, proportional transition rates from medication cycles to discontinuation periods and vice-versa. This group of assumptions/ restrictions conceptually simplify the recurrent multi-state structure back to a bidirectional-like structure which we refer to as “Emulated bidirectional” (Fig. 1G). Contrary to the previous structure (Fig. 1F) where the parameters of each transition were estimated separately, the joint parameter estimation among multiple transitions of the current structure (Fig. 1G) is computationally challenging and leads to excessive memory usage. Due to limitations in the maximum memory allocation (150GB), this structure was fit to a sub-sample of the study, including all the cases and randomly choosing 2 controls per case (1:2) while time-dependent effects of the case variable were allowed only for the transition from the starting state (start of follow-up) to the 1st medication cycle and death.

The estimated probabilities of the “Emulated bidirectional” structure (Fig. 1G) depicted in Fig. 7a, 7b, 7c and 7d can be interpreted in the same way as those of the recurrent multi-state structure without the restrictions (Fig. 1F), with the estimations of Fig. 6 being similar in shape and scale to those of Fig. 6. It can be observed that, contrary to Fig. 6, the estimated probabilities of Fig. 7 are very similar for cases and controls (with the exception of Fig. 7a). This is likely because we did not allow for time-dependent effects for the case variable among the transitions from medication cycles to discontinuation periods (and vice-versa), resulting in time-constant transition intensity rates ratios of cases versus controls that are close to 1 for those transitions (Figure A2b), leading to similar probability estimates of being in the same medication cycle among the two groups.

a) Total probability on a population level of being in a medication cycle since the start of follow-up, b) Probability of staying in current medication cycle since entering it. Total probability of being in a medication cycle since entering the 1st, 2nd, 3rd medication discontinuation period (d).

The advantages of emulating a bidirectional structure by applying restrictions to the initial recurrent structure is that we can tackle the issue of sparse data in high order states, gain precision in the estimations and conceptually simplify the multi-state structure back to a bidirectional-like structure, while allowing for different transition rates between each couple of medication cycle and discontinuation period states (and vice-versa). However, the “Emulated bidirectional” structure carries some of the drawbacks of the unrestricted recurrent multi-state model, making the assumption that individuals that enter the 6th medication cycle are characterized as chronic antidepressant medication users, while the initial bidirectional structure is free of this limitation. Additionally, due to its excessive memory usage, there are limitations as to how flexibly the baseline transition rates and the covariate effects can be modelled.