The Use of Real-World Data for Studies of Dynamic Disease Processes

2. Statistical models and methods for chronic disease processes

Some basics of multistate modeling. Methods of event history analysis⁴ are invaluable in the study of dynamic disease processes. In this article, we consider disease processes that can be conceptualized using multistate models wherein individuals may occupy different disease states over time; this can lead to insights concerning disease dynamics and related factors. Survival analysis can be viewed as a special case of multistate analysis where there are 2 states, “alive” and “dead”; however, even in survival analysis, > 2 states may be of interest when there are competing causes of death.⁵ Appreciation for the value of multistate models is growing in the public health arena; for example, there are applications directed at liver disease,⁶ end-stage renal disease,⁷ cancer outcomes following treatment,⁸ and dementia.⁹ In PsA, they have been used to model the cumulative number of damaged joints,¹⁰ the onset of axial disease,¹¹ or even progression of damage at the individual joint level.¹² Multistate models have also been used to study comorbidities of various diseases; for example, Perera et al¹³ studied changes in cognition in patients with systemic lupus erythematosus.

Multistate models are formulated by first specifying a set of states, 𝒮, which represent different stages of disease; for convenience, these are numbered 0, 1, 2, …, K. We will consider models for disease state as a function of time t since disease onset, with age of onset as a possible covariate of interest. We let Z(𝑡) denote the state occupied at time t ≥ 0, where Z(𝑡) takes a value in 𝒮 = {0, 1, …, K}. Multistate diagrams are used to depict a process, with arrows showing the pairs of states between which direct transitions are possible. Figure 1A shows a 4-state process, where state 0 may represent a healthy state, state 1 an initial diseased state, and state 2 a more advanced disease state; state 3 may represent death and is called an absorbing state. Such a model can, for example, be used to characterize the onset of psoriasis (entry to state 1) and the development of PsA among individuals with psoriasis (entry to state 2), while recognizing that individuals are at risk of death from each of states 0, 1, and 2.^14,15 Figure 1B is a simple progressive process that could be used to model the development of severely damaged joints as a function of the time t from the onset of arthritis; here, 𝑍(𝑡) = 𝑘 if an individual has accumulated 𝑘 severely damaged joints by time t. If one considers arthritis mutilans as having ≥ 5 severely damaged joints, we could set 𝐾 = 5 if we just want to model its occurrence. Although we could model the time to arthritis mutilans using a simple 2-state model (< 5 damaged joints and ≥ 5 damaged joints), the multistate approach here provides a more detailed and useful representation of the accumulation of joint damage. Note that there is an implicit death state in processes like Figure 1B that can be entered from any of the other states; when mortality rates are negligible, this is often omitted for simplicity.

• Download figure
• Open in new tab
• Download powerpoint

Figure 1.

Examples of multistate models, where the arrows depict the transitions that are possible at an instant in time. (A) A progressive multistate process with a terminal state. (B) A multistate representation of a counting process.

In multistate analyses, we formulate models for the probabilities of transitions between disease states as functions of time and past disease and covariate history. We let (𝑡) = {𝑍(𝑠), 0 < 𝑠 < 𝑡} denote the history of the disease process at time t, and V denote a vector of covariates (eg, demographic or genetic variables, age of disease onset). For convenience, we let 𝐻(𝑡) = {𝑍(𝑠), 0 < 𝑠 < 𝑡, 𝑉} denote the more complete history, including the covariate vector V. Intensity functions characterize the instantaneous risk of a transition (change in disease state) at any time t, given the history 𝐻(𝑡), and can be viewed as generalizations of hazard functions used in survival analysis.¹⁶ With multistate processes, however, intensity functions must be specified for all pairs of states for which instantaneous transitions are possible. The notion that the risk of events in the short-term future may be influenced by the past is natural in medicine, and intensity-based modeling has been used very effectively in rheumatology.^10,17,18

For a multistate process like the one in Figure 1A, the intensity for a 𝑘 → 𝑙 transition is the instantaneous probability that a 𝑘 → 𝑙 transition is made at time t, given state k is occupied just before time t (denoted as 𝑡⁻) and the previous process history. It is defined mathematically, as seen in Equation 1:

where 𝑡 + Δ𝑡⁻ means an instant before 𝑡 + Δ𝑡. The risk of transitions between disease states may depend on fixed covariates V (eg, sex, HLA markers), which are included in 𝐻(𝑡); the fixed covariates can explain heterogeneity in disease courses across individuals. In PsA, for example, women are known to experience a different disease course than men.¹⁹

Although Equation 1 may appear foreign, consider the case in which K = 1 and the state space is {0,1}, with only 0 → 1 transitions possible (ie, the 1 → 0 transition intensity is set to 0). In this case, 𝜆 (𝑡 ∣ 𝐻(𝑡)) reduces to the familiar hazard function for the time spent in state 0, which is commonly used in survival analysis. In Cox regression, for example, we specify multiplicative models where 𝜆₀₁(𝑡 ∣ 𝐻(𝑡)) = ℎ₀(𝑡)exp(𝑉′𝛽), with ℎ₀(𝑡) as a baseline hazard function and exp(𝑉′𝛽) a multiplicative factor reflecting the covariate effects. Specifically, exp(𝛽r) is the hazard ratio associated with a 1-unit difference in the rth covariate 𝑉r, with all other covariates held fixed.¹⁶

In practice, we formulate multistate models that express how the covariates and disease process history 𝐻(𝑡) influence change in disease states over time. For a pair of states k and l of the process in Figure 1, the k → l transition intensity would be of the form in Equation 2:

where 𝑉= (𝑉₁,…, 𝑉_p)′ is the set of covariates and 𝛽_kl = (𝛽_kl1,…, 𝛽_klp)′ is the corresponding set of regression coefficients. The function 𝜆_kl0(𝑡 ∣ (𝑡)) is called a “baseline” intensity, as it describes how the risk of a k → l transition depends on time 𝑡 and the disease history; note that the first subscript (k) in 𝜆_kl0(𝑡| (𝑡)) represents the state occupied at 𝑡⁻, the second subscript (l) represents the potential entry state under consideration, and the third subscript (0) in 𝜆_kl0(𝑡 ∣ (𝑡)) conveys that this is a baseline intensity function corresponding to an individual with 𝑉 = 0. As in Cox failure time models, exp(𝛽_klr) is interpretable as an instantaneous relative risk of a k → l transition for 2 individuals whose covariate value 𝑉_r differs by 1, while all other covariates are held fixed. Gladman and Farewell¹⁷ report such relative risks when studying the role of HLA markers on joint damage in PsA.

We previously provided details²⁰ on how to fit multistate models to data on a set of individuals and provided numerous examples, along with references to software. Cox model software can be used in special cases where exact times of transitions are always observable. However, state occupancy is often observed only at discrete times; for example, joint damage assessment requires radiographic examination that is done at periodic visits to a clinic. Markov models involve the simplifying assumption that 𝜆_kl0(𝑡 ∣ 𝑍(𝑡⁻) = 𝑘, (𝑡)) = 𝜆_kl0(𝑡), so that for the baseline intensity, the only relevant part of the disease history is time 𝑡 since onset and the state occupied at 𝑡⁻. This allows transition probabilities for successive intermittent observation times to be calculated without information on transitions between those times.²⁰ Most non-Markov models do not allow this to be done and are difficult to fit. For Markov models, the R package msm() enables model fitting and facilitates inference when observation is intermittent²¹; numerous illustrations are provided in Cook and Lawless.²⁰ Markov models are also easy to use for estimation of state occupancy probabilities as a function of time, which are useful for prediction. For example, estimates of 𝑃(𝑍(𝑡) = 𝑘 ∣ 𝑍(0) = 0, 𝑉), 𝑘 = 1,… can be used to predict the accumulation of joint damage t years after onset of PsA. For settings we consider in which processes are under intermittent observation, it is common and convenient to adopt parametric (including piecewise-constant) models for baseline intensities; the msm() package requires such a specification. If processes are subject only to right-censoring, semiparametric intensity-based models can be fitted²⁰; however, we do not consider this here.

Multistate models represent the dynamic aspects of a disease process in continuous time. Other time-varying factors such as biomarkers, visit times, and treatments associated with disease can also be included in models, as we discuss next. This is important for addressing challenges created by disease-related treatment decisions and observation times.

Joint models for markers and disease processes. For many diseases, there are dynamic serological biomarkers reflecting aspects of disease activity that are associated with disease progression. An example in arthritis is C-reactive protein (CRP), a systemic biomarker with high values occurring during periods of active inflammation. We let 𝑋(𝑡) denote the value of such a biomarker at time 𝑡 and let (𝑡) = {𝑋(𝑠), 0 < 𝑠 < 𝑡} denote the marker history. Jointly modeling 𝑍(𝑡) and 𝑋(𝑡) and can help in understanding the relationship between a dynamic marker and the disease process. We consider, for simplicity, a discrete marker, which could be obtained by categorizing continuous values according to whether they are low, moderate, or high, or simply elevated vs normal. In Section 5, we describe analyses where we dichotomize CRP as elevated or not.

In what follows, we consider a progressive disease process, in which the clinical states of interest are represented. Marker processes, for example, representing the level of disease activity or inflammation, are allowed to be relapsing and remitting, rather than progressive. Figure 2A is a multistate representation of a joint model for a binary marker {𝑋(𝑡), 0 < 𝑡}, taking the values of 0 (normal) or 1 (elevated), and a simple 2-state failure process {𝑍(𝑡), 0 < 𝑡}, where a 0 → 1 transition occurs upon failure (defined as the occurrence of some event of interest). Data on the marker are obtained only at discrete visit times, and we use gray arrows to indicate that we do not observe the precise times at which transitions are made between states 0 and 1, nor indeed the exact number of transitions. The times of some disease events (eg, failure) are observable; the black arrows for the 0 → 1 transitions of the 𝑍(𝑡) process convey this. An example is the time of occurrence of a cardiovascular event or the diagnosis of a comorbidity. In other cases, we may see that failure has occurred only at visit times. The effect of an elevated marker on the failure process can be ascertained by comparing the 0 → 1 intensities for the failure process when X(𝑡) = 1 vs when X(𝑡) = 0. Cox regression models are commonly used for this purpose and can be represented by the form of Equation 2 by adding 𝑋(𝑡⁻) to V. We remark that there are technical details concerning the rigorous mathematical treatment of such processes that we will not highlight here, as they have no effect on the methods we discuss. These details refer to the assumption that ≥ 2 simultaneous transitions are impossible, and to assumptions concerning whether transitions at a given time t are assumed to occur just before time t (that is, at time t⁻) or exactly at time t.

• Download figure
• Open in new tab
• Download powerpoint

Figure 2.

Examples of joint models for (A) a binary marker X(𝑡) and failure process Z(𝑡), and (B) a binary marker X(𝑡) and progressive disease process Z(𝑡) under intermittent observation. The dark arrows indicate transitions whose exact times are observable.

Figure 2B represents a more general progressive disease process {𝑍(𝑡), 0 < 𝑡} with state space {0,1,…,𝐾}. Such a process is depicted in isolation in Figure 1B, but in Figure 2B it is considered jointly with the marker process; a partial representation of this process is depicted for states 𝑘 − 1, 𝑘, and 𝑘 + 1. If 𝑍(𝑡) records the cumulative number of severely damaged joints at time t, for example, exact transition times are unknown, since this count is only observable upon clinical or radiographic examination; the gray arrows to the right convey this. The intensities in a parametric model can nevertheless be estimated, and we can assess the relative risk of disease progression for those with elevated vs normal markers (see our previous work²⁰).

As Figure 2 indicates, joint models are formed by defining states representing both the marker and disease process values, and we let 𝒵(𝑡) = (𝑍(𝑡), 𝑋(𝑡)) denote the state for the joint process. Thus, if {𝑍(𝑠), 0 < 𝑠} has 𝐾 + 1 states in Figure 1B, since the marker process can be in either of 2 states for every state of 𝑍(𝑡), the state space for the joint disease-marker process comprises 2(𝐾 + 1) states. Transitions are governed by the 2 types of intensities: one for the disease process and one for the marker process. We let ℋ(𝑡) = {𝒵(𝑠), 0 < 𝑠 < 𝑡, 𝑉} denote the history of the joint process; the transition intensities for the disease process {𝑍(𝑠), 0 < 𝑠} can then be defined as in Equation 3:

This allows the risk of transition between different disease states to depend on both current and past values of the disease and marker processes, as well as fixed covariates. Likewise, the intensities for the marker process are defined in Equation 4:

Taken together, Equation 3 and Equation 4 fully characterize the joint disease-marker process, such as the one represented in Figure 2B. The term “joint model” is used because the intensities for the disease process (Equation 3) and for the marker process (Equation 4) each depend on the histories of both processes. Specifically, this allows the risk of transition between disease states to depend on fixed covariates as well as past values of the disease and marker processes; likewise, the marker process transitions can depend on covariates and the history of the disease-marker process. Through the intensities, we can examine the effects of covariates V on both processes, calculate effects of markers on the disease process, and make predictions about the future disease course. We give specific examples of how intensity-based models for Equation 3 and Equation 4 can be set in the illustrative analyses of Section 5 (Equation 12 and Equation 13, respectively). We further expand this joint model to discuss modeling of treatment decisions in the next section.

Modeling treatment as a stochastic process. Consider patients attending the UTPAC and the prescription of biologic therapy. Individuals with an adverse risk profile or in an active disease phase may be more likely to be prescribed biologics than an individual whose disease is under control with milder therapies. The term “risk profile” refers to fixed personal attributes V that might be associated with a poorer prognosis; examples include early age of onset, family history, and HLA markers. An “active disease phase” is one in which pain is severe and inflammatory markers (eg, CRP) are elevated. Disease activity often has a relapsing-remitting nature, as reflected in the marker process {𝑋(𝑡), 0 < 𝑡} in Figure 2. If (𝑋(𝑡), 𝑉) are associated with treatment recommendations and decisions, and if they are not considered in analyses aiming to estimate a real-world effect of treatment, inferences are biased due to confounding by indication.²²

We consider the issues in a simplified setting to illustrate key points. We assume that biologics may be prescribed at any time after clinic entry and that once a prescription is made, an individual remains on some form of biologic therapy thereafter. If B is the time since disease onset of biologic prescription, we let 𝐵(𝑡) = 𝐼(𝐵 ≤ 𝑡) indicate treatment status at time t. We also let 𝑎₀ < 𝑎₁ < 𝑎₂ < ⋯ denote the times of clinic visits at which the disease state and marker measurements are made, where 𝑎₀ is the time at recruitment to the clinic. We let represent the number of clinical assessments to time t. We now define 𝒵(𝑡) = (𝐴(𝑡), 𝐵(𝑡), 𝑍(𝑡), 𝑋(𝑡)) as the state of the joint process that includes visits, biologic therapy, the disease process, and the marker process, and let ℋ(𝑡) = {𝐴(𝑠), 𝐵(𝑠), 𝑍(𝑠), 𝑋(𝑠), 𝑎₀ < 𝑠 < 𝑡, 𝑉} denote the complete process history to time t. Because the disease and marker processes are only observed intermittently, the observable history of the process is different; it is ℋ∘(𝑡) = {𝐴(𝑠), 𝐵(𝑠), 𝑎₀ ≤ 𝑠 < 𝑡, (𝑍(𝑎_j), 𝑋(𝑎_j)), 𝑗 = 0,…, 𝐴(𝑡⁻), 𝑉}. Note ℋ∘(𝑡) records the recruitment time, visit times, and the time biologics were prescribed if they were given over (0, 𝑡) along with the marker and disease states occupied at clinic visits.

The time 𝐵 that biologic therapy is prescribed is random and influenced by the disease process and marker values. This is modeled in Equation 5 using a corresponding intensity function:

Since this intensity accommodates a dependence on the history ℋ(𝑡), which includes the disease and marker process histories, it allows for the decision to prescribe therapy to be influenced by both features. As physicians aim to slow progression, reduce inflammatory markers, and treat based on clinical history, this aspect is appealing.

Figure 3 includes biologic treatment status, with the red arrows representing transitions that occur when biologic therapy is prescribed, for the joint model of Figure 2A. In Figure 3, since we are interested in the effect of biologics on failure, we do not consider cases where a biologic is prescribed after failure has occurred. Individuals may start biologic therapy at any time, and since the date that a prescription is filled is available, these transition times are known; the solid red lines convey this. If prescription of a biologic is driven only by information collected at the times of visits, one may assume that the prescription is made based only on observed data (ie, 𝜆^b(𝑡 ∣ ℋ(𝑡)) = 𝜆^b (𝑡 ∣ ℋ∘(𝑡))), but we do not need this assumption in joint models.²³

• Download figure
• Open in new tab
• Download powerpoint

Figure 3.

A joint multistate model for a binary biomarker, the prescription of biologic therapy, and a failure time. The dark arrows and the red arrows denote transitions whose exact times are observable.

Visit times may be scheduled according to a protocol, but, in reality, are often random and influenced by a person’s disease status. Equation 6 shows a model for the visit process expressed through an intensity:

6

This model allows for different individuals to have different propensities to attend the clinic according to covariates V, and also for visits to depend on the marker, disease, and treatment processes.

In the most general setting, the visit intensity at time t may be influenced by the disease, marker, and possibly treatment at 𝑡⁻, regardless of whether they are known. The formulation of a full joint model, and in particular, the Markov assumption, enables fitting of such a visit intensity. In much of what follows, we assume , which means that the timing of the jth visit can be influenced by the information observed up to the (𝑗 − 1)st visit, but not on any features of the disease process between visits. When the visit intensity satisfies this constraint, we refer to it as a conditionally independent visit process (CIVP). An expanded joint model can accommodate violations of this condition, and in such settings, classification of visits as either scheduled or disease-driven can be helpful.²⁴

This framework for visits allows for the construction of a joint multistate model for the disease, marker, treatment, and visit processes. Figure 4 depicts a generalization of the process in Figure 2A that incorporates the visit time process; transitions corresponding to clinic visits are represented by downward arrows. As before, the transition times corresponding to the black arrows can be observed precisely, and the transitions between the marker states (in gray) are not observable. Transitions depicted by the red arrows are observed, since the prescription date for biologics is known. Together, the intensities (Equation 3, Equation 4, Equation 5, and Equation 6) define the joint process in Figure 4.

• Download figure
• Open in new tab
• Download powerpoint

Figure 4.

An expanded joint multistate model for clinic visits, the prescription of biologic therapy, a failure time, and a binary biomarker. The dark arrows, red arrows, and green arrows denote transitions whose exact times are observable; green arrows represent new visits.

We note that for the Z and X processes, the model intensities at a given time 𝑡 depend on the state occupied at 𝑡, which can change between visits. Non-Markov models (eg, ones for which some intensities may depend on the time since state entry or other features of the process history) are difficult to fit in such settings,²⁰ and general software for dealing with them does not exist. Markov models can, however, be readily fitted based on the information observed at visit times using the msm() package. Model fitting is further simplified under the assumption of a CIVP.

The intensities in Equation 3, Equation 4, and Equation 5 define the joint model and based on these, a likelihood can be constructed for the observed data. Kalbfleisch and Lawless²⁵ describe how to do this and provide an algorithm for its optimization; the parameter values that optimize the likelihood are maximum likelihood estimates, which are consistent and normally distributed in large samples.²⁶ We discuss estimation further in Section 5 (Equation 15), with an explicit likelihood provided under the CIVP assumption.

Finally, the situation where 𝑍(𝑡) and 𝑋(𝑡) are measured at the same visit times is quite common, but if one of the disease or marker processes is not measured at certain visit times, the likelihood function must recognize this. This is a topic of current research, and we discuss it briefly in Section 6.

Models like those we have discussed are a plausible representation of many real-world processes involving the collection and study of data on disease. However, analysts often ignore the effects of markers and other variables on the prescription of treatment and the timing of observations. In the next section, we illustrate the effect of using simplified models that do not incorporate such potential confounding factors on the estimation of treatment effects.

3. Illustration of confounding by indication and model misspecification

In the previous section, we described a framework with which features of RWD can be emulated by considering dynamic disease and marker processes, a visit process, and treatment. We describe in an illustration in Section 5 how to fit models addressing these complexities, but in many applications, inappropriately simplified models are adopted. Here we consider the consequences of fitting a simple Cox regression model to examine the effect of a treatment on the time to a clinical event (failure) of interest. We consider both simple models examining only the effect of treatment as well as more elaborate ones that might be considered in recognition of the observational nature of the data and the risk of confounding by indication. All such models are misspecified since they do not recognize the complexities and the need for joint modeling.

We begin by fully characterizing the data generating process. We let V be a binary baseline covariate taking the value 1 with probability 0.5 and 0 otherwise. We consider event-free and biologic-naïve individuals as entering a clinic at age 𝑎₀ with a normal marker value. The changes in disease, marker, and treatment following recruitment are governed by the intensity functions that follow.

The intensity for the prescription of biologics (Equation 5) is specified as Equation 7:

7

So, the propensity to prescribe a biologic depends on the current status of the inflammatory marker (elevated vs normal). We note that 𝑅𝑅(𝐵|𝑋) = exp(𝛽₁) is the relative risk of biologics initiation at any given time when 𝑋(𝑡) = 1 (marker is elevated) vs when 𝑋(𝑡) = 0 (marker is normal). In what follows, we consider the influence of 𝑅𝑅(𝐵|𝑋) on the large sample estimate of the effect of biologics on failure risk obtained by fitting Cox regression models; in particular, the values of this parameter are represented on the horizontal axis for each panel of Figure 5.

• Download figure
• Open in new tab
• Download powerpoint

Figure 5.

Contour plots of the large sample percent relative bias in estimates of the effect of biologic therapy on failure risk from naïve Cox regression based on (Equation 10) to data generated under a conditionally independent visit process with (A) 𝜃₂ = 𝜃₃ = 0, (B) 𝜃₃ = 0, or (C) the full model in Equation 10, and (D) the full Cox model under a conditionally dependent visit process (non-CIVP); the true failure intensity is given in Equation 9 with 𝑒𝑥𝑝(𝜁₁) = 0.5, 𝑒𝑥𝑝(𝜁₂) = 2, and 𝑒𝑥𝑝(𝜁₃) = 1.5.

The marker alternates between 2 states and we let , k = 0,1, with Equation 8:

and

8

for the 1 → 0 and 0 → 1 transitions, respectively. The marker dynamics are influenced by V (if 𝜉₀₂ ≠ 0 or 𝜉₁₂ ≠ 0) and treatment (if 𝜉₀₁ ≠ 0 or 𝜉₁₁ ≠ 0). We set 𝜉₁₂ = log1.2 and 𝜉₀₂ = log0.8, so that more time is spent in the elevated marker state if 𝑉 = 1 compared to when 𝑉 = 0. Moreover, we let 𝜉₁₁ = − 𝜉₀₁ < 0 so that biologic treatment reduces the intensity for a transition to the elevated marker state and increases the rate of resolution if the marker is elevated; we let 𝑅𝑅(𝑋|𝐵) = exp(𝜉₁₁) and consider the influence of this on bias from naïve Cox regression analyses; this parameter is represented on the vertical axis of the panels in Figure 5.

We consider a failure process intensity (Equation 3) of the form in Equation 9:

9

where 𝜁₁ = log0.5, so treatment induces a 50% reduction in the intensity for failure given (𝑋(𝑡), 𝑉). We let 𝜁₂ = log2 characterize the effect of the marker on the failure process and let 𝜁₃ = log1.5, so that the intensity for failure is 50% higher for those with V = 1 compared to those with V = 0 given 𝐵(𝑡), 𝑋(𝑡).

We consider the process over the interval (0, 𝜏] with 𝜏 = 1 and assume 𝑎₀ = 0. We adopt time-homogeneous transition intensities and consider the 1 → 0 baseline intensity as 3 times larger than the 1 → 0 intensity for the marker process. We consider a probability of failure by the end of follow-up as 0.40, the probability of ultimately being prescribed biologics as 0.80, and specify an average of 8 visits by time 𝜏. For simplicity, for Equation 6, we focus on a visit process with 𝜆^𝑎(𝑡|ℋ(𝑡)) = 𝜆^𝑎, which means that visits are random in time but unrelated to the disease process or fixed covariates. A dependent visit process is obtained by setting 𝜆^𝑎(𝑡|ℋ(𝑡)) = 𝜆^𝑎exp(𝛾𝑋(𝑡)) with 𝛾 = log1.5, so individuals with an elevated marker have a higher visit intensity. The process exemplifies the complexities of RWD by defining a structure in which fixed covariates and a dynamic biomarker {𝑋(𝑡), 0 < 𝑡} affect risk of a failure time. The full process also incorporates the effect of fixed covariates and the dynamic biomarker on treatment decisions; the marker impact on the decision to prescribe is reflected by 𝑅𝑅(𝐵|𝑋). Once prescribed and taken, the treatment directly affects both the dynamic biomarker through 𝑅𝑅(𝑋|𝐵), and the disease process through the risk reduction 𝜁₁ = log0.50 in Equation 9.

Note that the full model also describes a setting in which time-dependent confounding²⁷ can arise from fitting simplified models because biologic-naïve individuals with 𝑋(𝑡) = 1 have a greater propensity to start biologics than those with 𝑋(𝑡) = 0. The marker {𝑋(𝑡), 0 < 𝑡} not only increases the likelihood of a biologic being prescribed but is also responsive to biologic therapy, in that individuals receiving biologics spend a reduced amount of time in the elevated marker state. Moreover, being in the elevated marker state increases risk of failure, so there is a consequential reduction in the overall risk of failure due to both the direct effect 𝜁₁ and the effect through the marker characterized by 𝑅𝑅(𝑋|𝐵) and 𝑅𝑅(Z|𝑋).

Here we consider the implications of adopting inappropriately simplistic Cox models in the face of such complexities, and we focus on the bias of treatment effect estimates. Our interest lies in illustrating the effect of confounding by indication²⁸ as well as that of model misspecification, since these 2 aspects are inextricably linked in the Cox regression model.²⁹ Using large sample theory on properties of estimators from misspecified Cox regression models,³⁰ we examine the percent relative bias of the effect of biologic therapy on the risk of failure from Cox models of the form of Equation 10:

A special case is Model A, in which 𝜃₂ = 0 and 𝜃₃ = 0, so that the only covariate is the time-dependent indicator of biologic therapy. Model B is a Cox model with 𝜃₂ = 0, so adjustment is made for the baseline covariate V only. Model C is the full Cox model in Equation 10 that also adjusts for 𝑋^∘(𝑡), the most recently recorded value of the marker 𝑋(𝑡) at 𝑎_A(t-). Note that although the true failure intensity (Equation 9) has a similar form to Model C, Model C conditions only on the observed (most recently recorded) version of the marker, rather than the true value, and so it is still misspecified.

Having described the complex interplay between the disease, marker, and treatment processes, we now consider the implications of fitting naïve Cox models on estimates of the effect of biologic therapy on failure. The percent relative bias of a large sample estimate for the effect of biologic therapy on failure risk from Cox models of the form (Equation 10) are evaluated in comparison to the true effect of biologics on the failure intensity given by 𝜁₁ in Equation 9; this represents the percentage of underestimation or overestimation of the real effect on the failure intensity. The bias will depend on several factors. We examine the bias as a function of the marker effect on the propensity to prescribe biologic therapy (𝑅𝑅(𝐵|𝑋)) and the effect of the biologic therapy on the marker process (𝑅𝑅(𝑋|𝐵)) when the association between the marker and the failure process is fixed at 𝑅𝑅(Z|𝑋) = 2; here the intensity for failure is doubled in individuals with an elevated (compared to normal) marker. Figures 5A-C give contours for the large sample percent relative bias for Models A-C as described here. Figure 5A shows that the effect of biologic therapy on the marker has the greatest impact, with the effect of the marker on the propensity to treat with biologics having a milder impact. Adjusting for the baseline covariate V has little effect on this bias (see Figure 5B vs Figure 5A). Adjusting for the most recently recorded marker value reduces the bias somewhat, but the trends remain. Adjusting for the true value of 𝑋(𝑡) and V would mitigate this bias, and this is achieved by fitting the joint model in Figure 3. Figure 5D shows how the large sample bias for Model C is further affected by a dependent visit process with 𝛾 = log1.5. In this framework, a joint model incorporating the visit process as depicted in Figure 4 would be required; we described this previously,²⁴ as did Lange et al.³¹

4. Assessing treatment effects with RWD

Contrasting RCTs and observational studies. The US Food and Drug Administration (FDA) defines RWD as “data relating to patient health status and/or the delivery of healthcare routinely collected from a variety of sources,” and RWE as “the clinical evidence about the usage and potential benefits or risks of a medical product derived from analysis of RWD.”³² The FDA documents its RWE program and lists ways in which RWE is currently being used. Much of this has been directed at pharmacovigilance studies to monitor drug use and side effects based on administrative data.³³ Efforts to use RWD for assessing effectiveness of therapeutic products have been more recent and often involve historical data. We take the term “evidence” to mean a body of results that can be reasonably interpreted as arising from a causal effect related to an underlying biological mechanism. We previously¹ noted differences between RCTs and observational studies when describing challenges in inferring the effects of treatment in real-world settings. We briefly review some of the issues, emphasizing the dynamic aspects of disease processes and treatment.

Clinical trials aim to examine the effect of experimental treatments on disease processes by focusing on patient-oriented clinical outcomes. A target population might be defined as the population of patients with a condition of interest (eg, PsA), but trial protocols usually define a narrower study population. Inclusion criteria might select individuals with active disease since they may provide evidence of a treatment effect.³⁴ For example, in a trial on biologic treatment, Mease et al³⁵ restricted attention to individuals with PsA who had ≥ 5 swollen and ≥ 5 tender joints along with a CRP level of ≥ 0.6 mg/dL. Once inclusion and exclusion criteria are set, the methods for identifying potential trial participants can be passive (eg, waiting to encounter suitable patients in a clinic) or proactive (eg, through solicitation from the population). For those in an RCT, dynamic marker and disease processes are in effect during the trial, but the protocol often defines a primary outcome that focuses on disease states or disease activity scores at a specified follow-up time 𝜏. Such outcomes can often be analyzed without a detailed dynamic model.

Relatively simple “marginal” measures of treatment effect are typically based on differences of simple averages or proportions. For example, Mease et al³⁵ use the change in disease scores during a trial period to define outcomes. Marginal effects are so named because they are defined by averaging (marginalizing) over known and unknown patient characteristics; randomization ensures that the distribution of these patient characteristics is the same in each treatment arm, so any differences in averages or proportions of outcomes seen in the respective arms can be attributed to the treatment they received. Such marginal effects tend to be less generalizable and transportable than effects that condition on observed individual-level covariates (see our previous paper¹ and Huang et al³⁶ for recent discussions). In the next subsection, we discuss marginal and conditional effects further, and in the following subsection, we discuss how both marginal and conditional measures of treatment can be based on dynamic joint models.

In contrast with RCTs, observational cohorts may be recruited from a target population through a poorly characterized process involving both active recruitment and self-selection by patients. Moreover, the standard of care is often not protocolized and may change over time, and a broader range of treatment options than in a typical RCT may be available. Patient assessment times during follow-up are often highly variable both within individuals over time and between individuals. Conceptualizing and defining an interpretable marginal treatment effect in settings involving a broad population, varied standard of care, treatment by indication, and irregular follow-up assessments is challenging. To mitigate the effect of confounding in assessing treatment effects, it is critical to understand the factors influencing decisions on when and how to treat patients. We address these issues through the models for disease processes and treatment discussed in Sections 2 and 3.

We note that the stochastic models we consider are consistent with broader views of causality and in particular, what is termed “dynamic causality.”⁴ This approach emphasizes plausible representations of process dynamics and understanding the effects of interventions on those dynamics. This can be done in the present context by considering the effect of treatment on the intensities for the disease and marker states. Dynamic causality is also related to the development of predictive models, which can be used to inform treatment decisions and policies. This framework is quite different than ones that define causal effects through potential outcomes³⁷; it takes the approach that a factor B (which may be time-varying) does not have a causal effect on an outcome Y if the inclusion of B in predictive models for Y (which will include other factors thought to be related to Y) does not increase predictive power.

Average, conditional, and personalized treatment effects. RCTs are primarily designed to estimate average treatment effects (ATEs), which—in the case of binary outcomes—are conventionally thought of as the difference in the marginal event probabilities for 2 treatment strategies. More formal causal reasoning expresses this effect as the mean difference in 2 counterfactual or potential outcomes for each patient: one applying if they receive the experimental treatment and the other if they receive standard care. Such a quantity is seemingly simple and easy to interpret, but factors that affect ATEs are often poorly understood for dynamic disease processes that evolve during the trial. Conditional treatment effects (CTEs) are defined by assessing differences between treatment groups among similar subjects. This can be achieved by examining treatment effects within regression models that adjust for subject-level covariates, or by forming subject subgroups or strata based on covariates. This is frequently done in RCTs under the heading of secondary subgroup analysis. It is also customary in multicenter trials to stratify on study centers to remove this component of variation. A key point in considering an ATE vs a CTE is that different types of individuals may respond differently to a given treatment; this is a fundamental concept of stratified medicine. Personalized medicine goes one step further and conditions on individual-level covariate values that may be unique to a given person; once again, this can be done using regression models. The use of very large numbers of covariates based on genomic data and other factors is a topic of much recent interest.³⁸

In the analysis of observational data, regression adjustment for fixed covariates V is crucial to mitigate the confounding effects of these covariates. We discuss ATEs and CTEs in observational settings in the next section, but we note here that CTEs that adjust for important covariates tend to be more transportable to other studies or populations. An ATE is typically less transportable because different studies or populations tend to have patients with different distributions of covariates.¹ In dynamic processes, which we address next, the same issues apply in further ways; for example, covariates may affect marker processes that are in turn associated with dynamic treatment processes. In this case, we need models like the ones described in this paper to make sense of treatment effects. Further, we need to understand the dynamics of treatment assignment in a real-world setting in order to consider alternative treatment plans or strategies, and how they can be compared.

Joint models for estimation of marginal and conditional effects in dynamic observational settings. A key issue for the specification and estimation of an ATE in real-world settings is the following: how should we define treatment “groups,” given that treatment changes can occur at arbitrary and disease-related times in the observational data? A common approach for the comparison of an experimental treatment and standard care^39,40 is to consider a sequence of “mini experiments,” defined using short time intervals (𝑡_j−1, 𝑡_j) and examine the set of individuals not on the experimental treatment at 𝑡_j−1. A comparison of outcomes over the interval is then made between those who initiate the experimental treatment during the interval and those who do not, with propensity score weighting or other adjustments for the dependence of treatment initiation on observed factors at time 𝑡_j−1. There are several problems with this approach, including that (1) short time intervals need to be used to minimize the effect of dynamic factors after time 𝑡_j−1 that affect treatment initiation by time 𝑡_j; (2) few persons will initiate treatment in a short time interval, so it is necessary to combine data across a set of intervals; (3) individuals can then appear in > 1 mini experiment, and can be in the standard care arm in one, and the treatment arm in another; and (4) it is cumbersome to combine the treatment effects across intervals so that they refer to an interpretable and relevant marginal outcome. Others, such as Gran et al,⁴¹ have considered the specification of a CTE using a similar approach; here we condition on confounding factors for each interval so that propensity score weighting is not required. However, the other problems remain. A final problem for such methods is that disease processes are under intermittent observation in most real-world settings. This means that the data needed to support such an analysis at a given time interval will typically be missing, making the procedures unworkable.

We propose here an approach that has been occasionally suggested for dynamic settings but not widely used to this point. The idea is to first specify a set of alternative treatment plans or policies; for convenience, we will use the term “plan” to refer to any of these. We then compare them in a hypothetical experiment in which the treatment plan is assigned randomly at a time 𝑡₀ for an individual, a schedule of assessments is set, and the individual is followed up to the time 𝜏 at which the outcome Y(𝜏) is observed. Treatment plans can be defined in various ways according to the actual real-world setting. For example, in the case of biologic treatment for PsA, 2 plans would be as follows: (1) remain off biologics over (𝑡₀, 𝜏); and (2) receive biologics over (𝑡₀, 𝜏). This corresponds to RCT designs used in the early assessment of new biolgics therapies. Regimens can, however, refer to ≥ 2 treatment plans and they can be of various types. They can be specified according to rules; for example, all individuals are on standard care at time 𝑡₀ and then (1) initiate biologics at the first time their marker 𝑋(𝑡) is observed in the elevated state, vs (2) initiate biologics if their marker has been in an elevated state for 2 consecutive visits. It is also possible to define treatment regimens by specifying that treatment intensities have specified forms (for examples, see Keiding et al⁴² and Gran et al⁴³). Once a set of alternative plans 𝑏 = 1,…𝑀 has been specified, we denote an individual’s assigned treatment plan by B^e, with superscript e indicating a hypothetical treatment regimen in which the treatment plan for an individual is randomly chosen from the alternative plans.

Then, we compute and compare outcome probabilities P^e(𝑌(𝜏) = 𝑦|B^e = 𝑏) for the different plans b = 1,2,…𝑀. These probabilities refer to the counterfactual outcome 𝑌(𝜏) that would occur if the individual was randomly assigned treatment plan b; we use the notation P^e along with B^e to refer to this hypothetical experiment. To estimate these probabilities, it is necessary to assume that the transition intensities for disease and marker states in the joint real-world process are the same as they would be in the hypothetical experiment. This is plausible when the intensities in the real-world fitted model condition on known confounders, and it defines the study population for the hypothetical experiment. To estimate an ATE based on 𝑌(𝜏), we can use what is sometimes called G-computation.³⁹ This requires that we first (1) fit our dynamic joint model to the RWD, conditioning on fixed covariates V and the initial disease and marker states at time 𝑡₀, and then (2) average over a distribution for V and the initial states using a specified distribution for them. This distribution could be the empirical distribution seen in the RWD, or another distribution specified to characterize the population in which the hypothetical experiment is conceptualized. We can similarly consider probabilities for 𝑌(𝜏) conditional on strata 𝑆 = 1,…,𝐾 based on subsets of V and initial states; comparisons of these probabilities for different treatment plans are a basis for CTEs. In this case, we would typically consider rather few strata in order to have a reliable estimate of the distribution of V and initial state, conditional on membership in stratum S = s.

We illustrate this approach in the next section.

5. Assessing the effect of biologics on persons with PsA

We illustrate the analysis of real-world disease process data and the approaches concerning treatment effects by considering biologic therapy for persons with PsA. We will use a comprehensive joint model for the biologics, disease, and marker processes, under the assumption that the visit process is ignorable, and we illustrate complexities of treatment assignment that are typical of clinical settings. We consider the specific goal of estimating the effect of biologic therapy on the risk of an increase of ≥ 2 severely damaged joints over a specified time interval. We consider time t as the time since disease diagnosis and the joint multistate model depicted in Figure 6 with states defined as follows. As earlier, we use 𝑋(𝑡) = 1 if the CRP marker is elevated and 𝑋(𝑡) = 0 otherwise. For the disease process, we first let N(𝑡) represent the number of clinically damaged joints at time t. Interest lies in the effect of biologics and so we restrict attention to the period in UTPAC from January 1, 2000, onward, since it was at this point that biologics began to be prescribed. We let v₀ be the date of the first clinic visit for a patient after this date, which occurs at time 𝑎₀ since disease onset; n₀ = N(𝑎₀) is the number of damaged joints at this point. A goal of biologic therapy is to reduce the progression of joint damage and so we define the disease process as {𝑍(𝑡), 0 < 𝑡}, let 𝑍(𝑡) = 0 if 𝑁(𝑡) = 𝑁(𝑎₀), 𝑍(𝑡) = 1 if 𝑁(𝑡) − 𝑁(𝑎₀) = 1, and 𝑍(𝑡)= 2 if 𝑁(𝑡) − 𝑁(𝑎₀) ≥ 2, for 𝑡 > 𝑎₀. This is the relevant process for estimating the effect of therapy on the risk of developing ≥ 2 additional damaged joints. Finally, we let B(𝑡) for those who initiated biologic therapy by time 𝑡 and B(𝑡) = 0 otherwise. Note that most individuals remain on some form of biologic therapy once they start such a therapy; we do not consider discontinuation or switching of biologics therapies, although one could do that if, for example, they were interested in comparing the effectiveness of different kinds of biologic therapy. This could be achieved by expanding the multistate diagram to allow for more states based on distinct therapies, but we do not pursue that here.

• Download figure
• Open in new tab
• Download powerpoint

Figure 6.

Joint model for 𝒵(𝑡) = (𝐵(𝑡), 𝑍(𝑡), 𝑋(𝑡)) for illustration. The red arrows represent treatment initiation; they are the only transitions whose exact times are observable.

Covariates may affect the disease process and treatment received and we let V represent covariates at 𝑎₀ including sex (male, female), BMI (calculated as weight in kilograms divided by height in meters squared), the Disease Activity Index for Psoriatic Arthritis (DAPSA; low, medium, high), the Psoriasis Area and Severity Index (PASI; mild, moderate, severe), pain as measured by the visual analog scale (VAS), and age at disease onset (years). Here DAPSA is coded as moderate if it was between 15 and 28 and high if it was > 28. PASI is coded as moderate if it was between 3 and 7 and high if it was > 7. We also include the number of damaged joints N(𝑎₀) as the extent of damage at 𝑎₀. We categorize this as follows at 𝑎₀ and subsequent times t: we let 𝑆 (𝑡) = 0 if 𝑁(𝑡) = 0, 1 if 𝑁(𝑡) ∈ {1,2,3,4}, 2 if 𝑁(𝑡) ∈ {5,6,7,8,9}, and 3 if 𝑁(𝑡) ≥ 10. Finally, we let 𝑆_k(𝑡) = 𝐼(𝑆 (𝑡) = 𝑘), where 𝑘 = 1,2,3 and define the vector S(𝑡) = (𝑆 ₁(𝑡), 𝑆 ₂(𝑡), 𝑆 ₃(𝑡))′.

There are a total of 921 individuals (56% male) contributing to this analysis, with a total of 15,800 visits and 10,862.21 person-years of follow-up. A total of 417 patients were ultimately prescribed biologics. At the baseline assessment, 301 (33%) had a moderate DAPSA and 364 (40%) had a high DAPSA, whereas there were 225 (24%) with a moderate PASI and 227 (25%) with a high PASI. The mean baseline VAS was 4.3 (SD 2.1) and the mean BMI was 29 (SD 6.1).

It is important to accommodate calendar time trends in the propensity to prescribe biologic therapy since its availability and use has steadily increased since 2000. To address this, we let 𝒞₀ < 𝒞₁ < ⋯ < 𝒞₅ denote dates (in dd/mm/yyyy format) separated by 5-year intervals with 𝒞₀ = 01/01/2000, 𝒞₁ = 01/01/2005, …, 𝒞₅ = 01/01/2025. Then, let 𝜏_j = max(𝒞_j − 𝐷, 0), 𝑗 = 0,1,…,5 be the disease durations at the respective dates, with 𝜏_j = 0 if disease onset was after 𝒞_j, j = 0,1, …,5. We then consider a defined time-dependent covariate 𝐶(𝑡) = 𝑘 if 𝑡∈ (𝜏_k, 𝜏_{k+ 1}), 𝑘 = 1,…,4, which helps distinguish the propensity to prescribe biologics over the 5-year windows of calendar time. Then, letting 𝐶_k(𝑡) = 𝐼(𝐶(𝑡) = 𝑘), 𝑘 = 1,…,4, we define 𝐶(𝑡) = (𝐶₁(𝑡),…, 𝐶₅(𝑡))′. The Lexis diagram in Figure 7 shows how the vector 𝐶(𝑡) distinguishing calendar time changes value in terms of the time since disease onset.

• Download figure
• Open in new tab
• Download powerpoint

Figure 7.

Lexis diagram illustrating the different time scales for intensity-based analyses, including the way calendar time trends in the tendency to prescribe biologics can be incorporated on the disease duration time scale.

With a right-censoring time R, the data for a particular individual can be written as 𝒟 = {𝐵(𝑠), 𝑎₀ ≤ 𝑠 < 𝑅, (𝑍(𝑎_r), 𝑋(𝑎_r)), 𝑟 = 0,1,…, 𝐴(𝑅), V}. With a set of m independent individuals, we have individual-level data 𝒟 = {𝒟_i, 𝑖 = 1,…,𝑚}, which we can use to fit the models described next.

To model the joint process 𝒵(𝑡) = (𝐵(𝑡), 𝑍(𝑡), 𝑋(𝑡)) in Figure 6, we specify the particular forms of the transition intensities. The individual-level models implicitly condition on 𝐵₁(𝑎_i0) = 0, 𝑋_i(𝑎_i0), 𝑁(𝑎_i0), V_i, where V_i represents covariates that may include relevant information on the disease history up to to 𝑎_i0.

We consider the intensity for biologics initiation as Equation 11:

11

which allows the propensity to treat a patient to depend on their inflammatory marker, fixed covariates, calendar time, and the extent of their joint damage. The k → k + l intensity (k = 0,1) for the joint damage (disease) process is modeled as Equation 12:

12

so the risk of progression is influenced by the marker status, biologic therapy, fixed covariates, and the cumulative burden of joint damage. The intensities for the marker process are modeled as Equation 13:

13

These are affected by biologic therapy, covariates, and the cumulative number of damaged joints. Finally, all the baseline intensities are assumed to be constant over time; this could be relaxed.

Under a CIVP assumption, the likelihood for a sample of n individuals is shown in Equation 14:

14

where is the observed history of the {𝒵_i(𝑠), 0 < 𝑠} process and the subscript 𝒵 for 𝑃_𝒵 means that this probability is computed based only on the intensities of the {𝒵(𝑠), 0 < 𝑠} process given earlier. For the Markov models, we consider this can be written simply as Equation 15:

15

(see Kalbfleisch and Lawless²⁵). This likelihood can be seen to involve only the observed data, but specification of the intensities (Equation 11, Equation 12, and Equation 13) and in particular, the Markov assumption they represent, means that the processes may change in continuous time between visits. The R function msm() by Jackson²¹ can be used to optimize this likelihood and obtain standard errors; we use it in what follows. If there is concern that the visit process may be related to the disease process, then a modified likelihood described by Lange et al³¹ and in our previous paper²⁴ must be used, but specialized code is required for analysis.

The results of fitting this full joint process are displayed in the Table. The first set of results is for the models for disease progression and for initiation of biologics. Here we see that the risk of damage progression is higher for patients with elevated CRP, higher DAPSA, and more damaged joints. We also see that older patients have a higher risk of damage progression. The intensity for disease progression is significantly reduced by the prescription of biologic therapy with a relative risk of 0.50 (95% CI 0.42-0.60; P < 0.01). There is significant increase in the intensity for the prescription of biologic therapy for those with an elevated CRP, who are male, and with higher DAPSA, high PASI, and higher pain scores on the VAS. We also see that, among patients with similar presentations, there is a reduced tendency for older patients to be prescribed biologic therapy. There is a significant effect of calendar time on the prescription of biologic therapy, with a significant increase in the 2006-2010 period compared to the 2000-2005 period, which persisted over the 2011-2015 and 2016-2020 periods. This can be explained by the increased availability of biologics and changes in clinical practice. The parameter estimates for the intensities of the CRP marker process are reported in the bottom half of the Table. There we can see that biologic therapy significantly reduces the risk that an individual with a normal CRP level transitions to an elevated value. Those with a higher PASI score, higher BMI, higher VAS, and many damaged joints have a higher risk of the CRP becoming elevated from a normal state; older patients presenting with the same condition have a reduced intensity for an elevated CRP. At times when CRP values are elevated, those with higher DAPSA and higher BMI have lower intensities for the transition to normal CRP values; those with a moderate number (5-9) of damaged joints have a reduced intensity for the normalization of CRP values compared to those with 0 damaged joints. The regression coefficient for biologic therapy surprisingly suggests those on biologics will have a reduced intensity for the normalization of the CRP marker; this may arise from residual effects of confounding by indication and the fact that the marker is dichotomized. In ongoing work, we are exploring the possibility of defining more marker states.

The insights from this fitted model are of great help in understanding the disease process. It is clear that there is a complex interplay between the inflammatory marker, the prescription of biologics, and the joint damage process. The effects of biologics in these models are expressed as multiplicative effects on intensity functions. These reflect an instantaneous effect on the risk of events, conditional on the process history, and so we have what is referred to as a dynamic causal interpretation. That is, if the intensity-based models are an accurate reflection of the process at an instant in time t, an individual on biologics is expected to have a different future course than an individual who is untreated; in the short term, this difference arises because the intensities governing change differ between the treated and untreated individuals. There is no conceptualization of potential outcomes here, so the coefficients do not correspond to a causal effect in this framework. For an excellent discussion of the different schools of causal reasoning see Aalen and Frigessi.⁴⁴ For estimation of an ATE of biologic therapy, which might be of interest when examining real-world effects to compare with effects seen in RCTs, we use a G-computation approach,⁴⁵ described in the preceding section. This involves using the fitted intensity functions under 2 hypothetical experimental treatment plans for each patient in the sample, and we describe this here. CTEs can also be considered, but we will not address them here.

The 2 treatment plans we consider are ones where an individual is either receiving biologics or not receiving biologics, respectively, for the entire follow-up period (𝑎₀, τ). For the former, we can calculate the conditional probability of an increase of ≥ 2 damaged joints, given baseline covariates, for each individual i = 1, …, m in the UTPAC data, using the estimated disease damage and marker intensities on the B(𝑡) = 1 side of Figure 6. We can then calculate the average effect by averaging over the m individuals in the cohort; this effectively uses the empirical distribution of the covariates in these individuals. Similarly, we can calculate the conditional and average probabilities for the regimen, where B(𝑡) = 0 for all persons, using the estimated intensities in the B(𝑡) = 0 side of Figure 6. We can then compare the estimated average probabilities of an increase of ≥ 2 damaged joints for each of the 2 treatment regimens, and define an ATE as the difference or ratio of the 2 probabilities.

In this hypothetical experiment, the probability of a counterfactual outcome when a person is given treatment b = 0 or b = 1 is, in effect, calculated for each individual. This can be computationally demanding, although methods and software for doing this exist.²⁰ Another approach is to use simulation; in this case, we simulate from the fitted joint multistate process, using the B(𝑡) = 0 and B(𝑡) = 1 sides of Figure 6, respectively. Specifically, for the case where B_i(𝑎₀) = 1, we simulate the joint process in Equation 16:

16

where the superscript reflects the fact that the path is simulated under this hypothetical treatment plan. For the second plan, where B_i(𝑎_i0) = 0, we simulate the joint process conditional on ℋ_i(𝑎_i0) = 1, which includes B_i(𝑎_i0) = 0 while setting the intensities for biologics initiation to zero; that is, we stay on the B(𝑡) = 0 side of Figure 6.

The “control” plan just described assumes that an individual cannot initiate biologics during the follow-up period; this may be viewed as representing standard care in the setting where biologics are not available. In some settings, especially when 𝜏 – 𝑎₀ is fairly large, it may be preferred to consider a trial that allows switches to biologic therapy during follow-up. In this case, we need to specify biologic initiation intensities that allow this; the default option is to use the estimated intensities from the real-world UTPAC data, but other choices could be made for the hypothetical experiment. Individual-level data for the control regimen would in this case come from the process in Equation 17:

17

We would use an intention-to-treat comparison here, defining the 2 treatment regimens according to individuals’ initial biologics status at 𝑎_i0. For the illustration that follows, we assume the control regimen does not allow a switch to biologics treatment during follow-up. We also employ simulation to obtain ATEs, as we describe below.

The UTPAC data used to fit the full joint process are based on intermittent observation of the disease and marker processes under the assumption of a CIVP.²⁴ However, when simulating data under the 2 regimens, we can generate the joint process in continuous time. This produces for each individual the binary counterfactual outcomes and which indicate whether or not there has been an increase of ≥ 2 damaged joints by time units after the initial visit. A simple estimate of a marginal effect (or ATE) of biologics is then , where , j = 0,1. We consider 𝜏_o = 2 years here as the value of prime interest but consider a range of values in the application that follows.

Note that this ATE corresponds to a population average effect in which the sample of individuals in the observational study is viewed as representative of the target population. One may alternatively average over a subsample of the UTPAC sample by applying inclusion criteria based on the disease state, marker state, and covariates; this may be done by confining attention to individuals with histories at their respective 𝑎_i0, satisfying certain criteria. For example, one may examine the effect of giving biologics earlier during the disease course by restricting attention to individuals with 𝑎_i0 less than some threshold, or examine effects in patients with more extensive damage at time 𝑎_i0.

For the illustration here, we used broad eligibility criteria and included all patients in UTPAC who are biologic-naïve at their first visit. The estimated marginal probabilities of developing ≥ 2 additional damaged joints are given as a function of τ in Figure 8A, along with symmetric pointwise 95% CI obtained by estimating the standard error using the parametric bootstrap with 100 samples. The 2 treatment plans are (1) where patients remain biologic-naïve over, (𝑎_i0, 𝑎_i0 + τ) and (2) when they are prescribed biologics at 𝑎_i0 and stay on biologics. The risk of developing ≥ 2 additional damaged joints is lower for the plan in which biologic therapy is given at 𝑎_i0 (at 𝑎_i0 + τ_o the risk is 0.025 [95% CI 0.018-0.031]) compared to the plan where it is not prescribed (here the risk is 0.08 [95% CI 0.07-0.09]). In Figure 8B, the difference in these estimated probabilities is plotted as a function of τ, where at 2 years there is a −0.05 (95% CI −0.06 to −0.04) absolute reduction in the risk that ≥ 2 additional joints become damaged corresponding to the use of biologic therapy.

• Download figure
• Open in new tab
• Download powerpoint

Figure 8.

(A) The cumulative risk of ≥ 2 additional damaged joints in untreated and biologic-treated patients, and (B) the average absolute reduction in risk of ≥ 2 additional damaged joints from biologic therapy. Pointwise 95% CIs based on the parametric bootstrap are also shown in light lines for both panels.

We reiterate the strong implicit assumption, mentioned in Section 4 on joint models for estimation of marginal and conditional effects, that the intensities estimated using the RWD apply to the hypothetical experiment (or trial) considered. This supposes (1) that all important covariates affecting all transition intensities have been identified, (2) that their effects have been correctly modeled, and more subtly, (3) that the intensities estimated in the UTPAC “real-world” are transportable to the hypothetical world in which treatment is assigned at random at the initial visit. A salient point is that a disease process is likely to differ in a controlled experimental trial and in a real-world observational setting where there is a less rigid protocol governing standard of care and compliance with the assigned treatment. This can be addressed to some extent by carefully considering these issues in the real-world setting and how the process intensities might differ from those in an RCT with the same factors observed.