Uncovering the heterogeneity of disease impact in axial spondyloarthritis: bivariate trajectories of disease activity and quality of life.

OBJECTIVE: The goal of managing axial spondyloarthritis (axSpA) is to improve and maintain patients' health-related quality of life (HRQoL), mainly through targeting towards low disease activity. Here, we aim to gain insight into the joint evolution of HRQoL and disease activity by identifying and characterising latent subgroups of patients with longstanding disease displaying similar trajectories throughout 8  years of follow-up.
METHODS: Data from Outcome in Ankylosing Spondylitis (AS) International Study (n=161) and Groningen Leeuwarden AS cohort (n=264) were used. Biennially, HRQoL was assessed by AS Quality of Life (ASQoL) and disease activity by AS Disease Activity Score-C reactive protein (ASDAS-CRP). Bivariate trajectories of these outcomes were estimated by group-based trajectory modelling. Next, trajectories were profiled by comparing the latent groups with respect to baseline factors using analysis of variance and χ² test.
RESULTS: Five bivariate trajectories were distinguished, in which ASQoL and ASDAS-CRP were tightly linked: (t1) low impact of disease; (t2) moderate impact; (t3) high impact with major improvement; (t4) high impact with some improvement; (t5) very high impact. Profiling revealed, for example, that (t1) was characterised by male gender and Human Leucocyte Antigen B27 positivity; (t3) by younger age, shorter symptom duration and biological intake and (t5) by the highest proportion of females.
CONCLUSIONS: We identified five bivariate trajectories of HRQoL and disease activity demonstrating a clear mutual relationship. The profiles revealed that both individual-related and disease-related features define the type of disease course in respect to HRQoL and disease activity in axSpA. This may provide clinicians insight into the differences among patients and help in the management of the disease.

The goal of managing axial spondyloarthritis (axSpA) is to improve and maintain patients’ health-related quality of life (HRQoL), mainly through targeting towards low disease activity.

Trajectories with respect to disease activity have been identified in patients with early axSpA.

We identified five latent groups displaying similar trajectories of HRQoL and disease activity in patients with established axSpA.

The profiles revealed that both individual-related and disease-related features define the impact of axSpA.

The profiles may provide clinicians insight into the differences among patients and help in the management of axSpA.

Introduction

The clinical presentation of axial spondyloarthritis (axSpA) is heterogeneous and can include inflammatory back pain, arthritis, dactylitis and enthesitis. Furthermore, extra-articular manifestations commonly occur.1 As a consequence, patients might experience pain, fatigue, impaired physical function and restricted participation in social roles.2 This results in reduced health-related quality of life (HRQoL) and lower life-satisfaction when compared with the general population.3 In addition, restrictions in work participation add to the socioeconomic impact of the disease for society.4 5 In clinical practice, improving and maintaining patients’ HRQoL is the ultimate goal of management of axSpA.6 Optimisation of HRQoL can be achieved through management of disease activity by means of (pharmacological) therapy, for which a treat-to-target strategy is increasingly advocated.6 7

While differences among individuals with axSpA exist with respect to the clinical presentation, the type and severity of disease characteristics within an individual may also change over time.8 9 Several studies tried to explain variability in HRQoL and disease activity in axSpA using classic regression analyses. Variability in HRQoL has been related to pain, stiffness and fatigue and also to personal contextual factors such as education, helplessness or comorbidities.7 10–12 Variability in disease activity has been related to age, Human Leucocyte Antigen B27 (HLAB-27) positivity, smoking and pharmacological treatment among which NSAIDs and tumour necrosis factor α (TNFα)-inhibitors.13–16 However, none of these studies investigated the presence of possible heterogeneity (ie, subgroups) within their population that is not necessarily attributable to known covariates. To gain more insight into this heterogeneity, a change of perspective might be of interest. Instead of addressing and explaining the variability of axSpA outcomes by the usual variables-centred approaches (like regression models), it is possible to take the variability itself as a starting point to cluster patients according to their clinical presentation over time. Multiple, not directly observable subgroups of patients are likely to emerge, each constituting a different symptomatic course of disease. These latent subgroups are represented by longitudinal trajectories, which capture interindividual differences in intraindividual changes of selected outcome(s). The heterogeneity of clinical presentation in axSpA warrants the application of such person-centred approach, mainly because substantial variability is expected in the (course of) impact of axSpA and effect of treatment. In addition, revealed subgroups can be described by means of patients’ characteristics resulting in risk-stratification that may more closely meet the demands of daily clinical practice.

To gain such understanding, Moltó et al evaluated in a French cohort of patients with early axSpA the presence of potential latent subgroups with respect to disease activity during a 3-year follow-up period.17 Five disease activity trajectories were observed: (1) persistent moderate disease activity; (2) persistent inactive disease; (3) changing disease activity; (4) persistent high disease activity; (5) persistent very high disease activity. Interestingly, in their study, a relation between the disease activity trajectories and the level of HRQoL at baseline was established. Both HRQoL and disease activity are factors that reflect the impact of the disease in a person. As management of disease activity should also be mirrored in gains in HRQoL, it is important to understand their joint evolution, that is, how temporal patterns of codependencies of the two outcomes unfold. The role of treatment herein, which is highly effective in controlling disease activity, would be of additional clinical relevance.

In the present study, we aimed to explore the heterogeneity of the impact of axSpA by identifying and characterising latent subgroups of patients with similar trajectories of HRQoL and disease activity in two well-phenotyped cohorts of patients with longstanding disease, who were followed biennially up to 8 years.

Patients and methods

Patients

Data from two prospective, multicentre, longitudinal observational cohort studies, Outcome in Ankylosing Spondylitis (AS) International Study (OASIS) and Groningen Leeuwarden AS (GLAS) cohort, were used for the present study.18 19 OASIS started in 1996 and was conducted at several secondary and tertiary referral centres in the Netherlands, Belgium and France. At baseline, 217 patients were included, all fulfilling the modified New York criteria (mNYC).20 In order to increase the sample size and create a database most related to reality with all types of disease activity and treatments, including anti-TNF-α treatment, which was not yet available at the start of OASIS, the sample was enriched with patients from the GLAS cohort. GLAS started in 2004 and was conducted in the north of the Netherlands. From this cohort, 266 patients who started anti-TNF-α treatment because of active disease between 2004 and 2012 were included. Patients fulfilled either the mNYC or the imaging arm of the Assessment in SpondyloArthritis International Society criteria for axSpA.20 21 Patients from both cohorts were followed, according to a fixed protocol at regular intervals and follow-up continued also in patients stopping/switching treatment. All patients provided written informed consent according to the Declaration of Helsinki.

Outcome variables

HRQoL was assessed at each visit by the AS Quality of Life (ASQoL) questionnaire.22 The ASQoL is a need-based HRQoL questionnaire consisting of 18 impairments and limitations typical for axSpA. Score ranges from 0 to 18, higher scores imply worse HRQoL. A threshold of 8 was used to define a patient acceptable symptom state.23

Disease activity was evaluated at each visit by several constructs that allowed to calculate the AS Disease Activity Score with C-Reactive Protein (ASDAS-CRP).24 ASDAS-CRP cut-offs were used to describe disease activity states: ‘inactive disease’ (<1.3), ‘low disease activity’ (≥1.3 and<2.1), ‘high disease activity’ (≥2.1 and≤3.5) and ‘very high disease activity’ (>3.5). A change of ≥1.1 in ASDAS-CRP was defined as ‘clinically important improvement’ and ≥2.0 as ‘major improvement’.25

The following baseline parameters were available to further characterise the latent subgroups: age, gender, smoking status (yes/no), body mass index, symptom duration, CRP-level, HLA-B27 (positive/negative), presence of (bridging) syndesmophytes (yes/no) on spinal radiographic imaging and start of anti-TNF-α treatment (yes/no). Further, Bath AS Disease Activity Index (BASDAI),26 Bath AS Functional Index (BASFI)27 and both patient and physician Global Disease Activity were assessed.

Statistical analysis

The analyses proceeded in two phases: (1) model selection consisting of identification of the latent trajectory groups and their functional forms and (2) characterisation of the unveiled groups. Patients without any ASQoL measurement were excluded from the analyses.

To identify clinically distinct trajectory classes of HRQoL and disease activity, group-based multitrajectory modelling (GBTM) was used, which is a form of latent class analysis.28 29 A major distinction between the multivariate GBTM, compared with its univariate version, is that in the former each identified trajectory is a representation of the temporal course of several outcomes jointly (herein specifically HRQoL and disease activity). Thus, the bivariate-trajectories in this study were estimated simultaneously, not individually, one at a time. As a result, each identified latent class captured individual patterns of changes of HRQoL and disease activity and the dynamics of their temporal codependencies. With GBTM, no prior defined classes are presumed, and therefore, heterogeneity in a population that is not attributed to any known factor can be uncovered. Models with increasing numbers of classes were run (from 1 up to 7) to determine the number of latent classes (class-enumeration) with SAS procedure PROC TRAJ.28 We resorted to the Fit-Criteria Assessment Plot (FCAP) R-code to assist the process of class-enumeration.29 FCAPs are graphical displays providing eight goodness-of-fit and model-adequacy criteria (Bayesian Information Criterion, Akaike’s Information Criterion, likelihood of assignment, average posterior probability of assignment, mismatch between estimated and assigned group probabilities, SD of group membership probabilities, odds of correct classification, proportion of individuals estimated to be assigned to the smallest group). After settling the number of latent classes, the polynomial order of the trajectories determining their final shape was selected, based on statistical significance (see online supplementary material for a detailed description of the model fitting procedure). GBTM estimates each participant’s probabilities for trajectories’ membership posthoc. A patient is then classified into the latent class for which (s)he shows the highest posterior probability of assignment.

Subsequent to model-selection, extracted HRQoL-disease activity trajectories were compared with respect to several baseline demographic and clinical variables with analysis of variance (posthoc Bonferroni correction) and χ² test, where appropriate.

Intermittent missingness (ie, a single variable or an entire visit missing, followed by complete assessment) was assumed to be missing at random. Loss to follow-up, that is, non-random attrition occurred. To test for differential drop-out, we computed for all classes the proportion of individuals with two or more missing values in the final and immediate preceding time points. Their frequencies were then compared.

Analyses were performed using SAS V.9.4, SPSS V.21 and R V.3.2.0. The significance level was set at 0.05.

Results

From the 217 patients included in the OASIS cohort, we removed the first 4 years of follow-up, because the ASQoL measurement was not available at that time. One patient was excluded from all (previous and current) analyses from OASIS, because of inconsistencies in the data. Another 56 patients had no ASQoL measurements after year 4. In total, 161 (74%) of the patients from the OASIS cohort were included in the current study. These patients did not differ in baseline characteristics from the entire cohort. From the GLAS cohort, 2 patients were excluded because of unavailability of ASQoL measurements, and 264 patients (99%) were included.

Ultimately, the analyses were conducted on a multicohort of 425 unique patients, who were followed biennially up to four time points after inclusion in the present study. Sometimes, individual patients missed measurement points. At time point 0, data from 411 patients were available, after 2 years from 363 patients, after 4 years from 290 patients, after 6 years from 238 patients, and after 8 years from 158 patients.

Identifying latent HRQoL-disease activity trajectory groups

Baseline characteristics of the overall patient population and separately for each cohort are provided in table 1. Patients included in OASIS were characterised by older age, longer symptom duration and higher proportion of (bridging) syndesmophytes in comparison to the GLAS population. Due to the inclusion criteria of GLAS, these patients were characterised by higher CRP, BASDAI and proportion prescribed TNFα-inhibitors.

Table 1

Baseline characteristics of the patients

Characteristics	OASIS(n=161)*	GLAS(n=264)*	P values	Overall(n=425)*
Age (years)	49.1 (11.7)	42.7 (11.9)	<0.001	45.2 (12.2)
Male sex	112 (69.6%)	179 (67.8%)	0.704	291 (68.5%)
Smoking—yes	41/116 (35.3%)	88/229 (38.4%)	0.576	129/345 (37.4%)
BMI	26.5 (4.9)	26.6 (4.3)	0.817	26.5 (4.5)
Symptom duration	26.0 (11.5)	17.2 (11.3)	<0.001	20.6 (12.2)
CRP (mg/L)	12.3 (17.0)	18.0 (18.3)	0.002	16.0 (18.1)
HLA-B27-positive	129/156 (82.7%)	204/257 (79.4%)	0.409	333/413 (80.6%)
Presence of (bridging) syndesmophytes	74/109 (67.9%)	107/197 (54.3%)	0.021	181/306 (59.2%)
Start anti-TNF-α treatment	1/154 (0.6%)	256/262 (97.7%)	<0.001	257/416 (61.8%)
BASDAI (0–10	3.3 (2.3)	6.1 (1.7)	<0.001	5.1 (2.3)
BASFI (0–10	3.6 (2.7)	5.6 (2.2)	<0.001	4.9 (2.6)
Patient GDA (0–10)	3.4 (2.7)	6.7 (1.9)	<0.001	5.5 (2.7)
Physician GDA (0–10)	2.0 (1.6)	4.7 (2.2)	<0.001	3.8 (2.4)
ASDAS-CRP	2.5 (1.1)	3.8 (0.8)	<0.001	3.3 (1.1)
ASQoL (0–18)	5.8 (4.7)	10.0 (4.1)	<0.001	8.7 (4.7)
Attrition ASQoL (≥2 missing final measures)	42 (26.1%)	57 (21.6%)	0.288	99 (23.3%)
Attrition ASDAS-CRP (≥2 missing final measures)	62 (28.5%)	56 (21.2%)	<0.001	118 (27.8%)

Values represent mean (SD) or number (%).

*Number of patients per group, unless otherwise stated due to missing values.

anti-TNF-α treatment, antitumour necrosis factor α treatment; ASDAS-CRP, Ankylosing Spondylitis Disease Activity Score—C reactive protein; ASQoL, Ankylosing Spondylitis Quality of Life; BASDAI, Bath Ankylosing Spondylitis Disease Activity Index; BASFI, Bath Ankylosing Spondylitis Functional Index; BMI, body mass index; CRP, C reactive protein; GDA, Global Disease Activity Score; GLAS, Groningen Leeuwarden Ankylosing Spondylitis; HLA-B27, Human Leucocyte Antigen B27; OASIS, Outcome in Ankylosing Spondylitis International Study.

According to model-adequacy criteria, the five-class solution provided the best fit to the data (online supplementary figure 1 and online supplementary table 1). The spaghetti plots (figure 1) display individual-specific lines within each of these five latent classes. Each extracted latent class is assigned one colour and represented by two developmental trajectories: an ASQoL (HRQoL)-CRP and an ASDAS-CRP (disease activity) trajectory. Therefore, the developmental paths displayed convey individual temporal progressions of HRQoL and disease activity and their coevolution. For instance, subjects assigned to class 1 (green) showed a gradual improvement in disease activity with concurrently a good and stable temporal unfolding of HRQoL. Similarly, those patients assigned to trajectory 3 (black) showed a major decrease in disease activity and had at the same time a substantial improvement of HRQoL, already after baseline. By analogy, the same interpretation of temporal covariation between disease activity and HRQoL extends to the other developmental classes. Since the trajectories illustrate concomitant patterns of change and/or stability of HRQoL and disease activity over time, the distinct trajectories are denoted: (t1) Low impact of axSpA on a patient: stable low and acceptable ASQoL (values ranged from 0.8 to 1.6 and stable ASDAS-CRP inactive to low disease activity (range 1.2 to 1.7); (t2) Moderate impact: stable moderate but acceptable ASQoL (range 4.5–5.6 and stable ASDAS-CRP high disease activity (range 2.3–2.8); (t3) High impact with major improvement: initial unacceptable ASQoL (point estimate 8.9) and ASDAS-CRP very high disease activity (point estimate 3.9), but major improvement in both (range 1.6–2.2 for ASQoL and 1.6–1.9 for ASDAS-CRP); (t4) High impact with some improvement: a minor improvement from unacceptable average ASQoL (point estimate 11.0) towards a lower but still unacceptable state (range 8.1–9.4), in combination with a clinically important improvement from ASDAS-CRP very high disease activity (point estimate 3.7) towards ASDAS-CRP high disease activity (range 2.1–2.6); (t5) Very high impact: persistently high and unacceptable ASQoL (range 12.6–14.4 with ASDAS-CRP very high disease activity (range 3.2–4.1). The estimated average trajectory lines of these five distinct developmental joint latent classes of ASQoL and ASDAS-CRP over the five time points are depicted in figure 1. For all trajectory-groups, HRQoL and disease activity mirrored to a great extent the temporal patterns of one another.

Figure 1

Joint trajectory latent classes of ASQoL–ASDAS-CRP. At the top: spaghetti plots of observed ASQoL and ASDAS-CRP values for each individual over a period of 8 years (five measurements), coloured according to their final latent class assignment. Note that each class is represented by both an ASQoL and an ASDAS-CRP trajectory. At the bottom: estimated average trajectories’ lines with corresponding 95% confidence bands and baseline prevalence. ASDAS-CRP, Ankylosing Spondylitis Disease Activity Score—C reactive protein; ASQoL, Ankylosing Spondylitis Quality of Life.

Characteristics of the trajectory groups

Baseline characteristics of patients in each latent class are shown in table 2. The trajectory of low impact of disease (t1) showed the largest proportion of males and HLA-B27 carriers, the lowest proportion of patients with (bridging) syndesmophytes and lowest proportion of patients prescribed anti-TNF-α treatment at baseline. In comparison to the low impact group, the moderate impact group (t2) is represented by a lower proportion of HLA-B27 carriers and a higher proportion of patients with (bridging) syndesmophytes. In the trajectory high impact with major improvement (t3), younger patients were found, having also the shortest symptom duration. Further, this group had the largest proportion of patients prescribed anti-TNF-α treatment. The trajectory high impact with some improvement (t4) had the highest proportion of (bridging) syndesmophytes as well as the longest symptom duration. Last, the trajectory very high impact of disease (t5) was represented by the highest proportion of females. In t4 and t5, similarly as in t3, a relative high proportion of patients was prescribed anti-TNF-α treatment at baseline, but without reaching the same effect in ASDAS-CRP and ASQoL.

Table 2

Baseline characteristics of patients within each HRQoL—disease activity trajectory

Characteristics	Trajectory 1Low impact(n=55)*	Trajectory 2Moderate impact(n=105)*	Trajectory 3High impact with major improvement(n=86)*	Trajectory 4High impact with some improvement(n=124)*	Trajectory 5Very high impact(n=55)*	P values
Patients from OASIS	40 (72.7%)	62 (59.0%)	4 (2.5%)	37 (29.8%)	18 (32.7%)	<0.001
Patients from GLAS	15 (27.3%)	43 (41.0%)	82 (95.3%)	87 (70.2%)	37 (67.3%)	<0.001
ASDAS-CRP	1.6 (0.7)	2.8 (0.8)	4.0 (0.7)	3.7 (0.8)	4.1 (0.7)	<0.001
ASQoL (0–18)	1.6 (2.0)	5.4 (2.7)	9.0 (3.5)	11.1 (2.8)	14.6 (2.4)	<0.001
Attrition ASQoL (≥2 missing final measures)	16 (29.1%)	17 (16.2%)	17 (19.8%)	32 (25.8%)	17 (30.9%)	0.146
Attrition ASDAS-CRP (≥2 missing final measures)	20 (36.4%)	26 (24.8%)	17 (19.8%)	36 (29.0%)	19 (34.5%)	0.158
Age (years)	43.9 (12.4)	46.0 (12.7)	40.9 (10.7)	47.9 (11.3)	45.1 (13.7)	0.001
Male sex	46 (83.6%)	71 (67.6%)	62 (72.1%)	83 (66.9%)	29 (52.7%)	0.012
Smoking—yes	16/42 (38.1%)	23/85 (27.1%)	26/74 (35.1%)	44/98 (44.9%)	20/46 (43.5%)	0.129
BMI	25.5 (3.7)	26.6 (4.9)	25.9 (3.9)	26.8 (4.3)	27.7 (5.4)	0.210
Symptom duration (years)	20.3 (10.9)	21.4 (12.1)	17.0 (10.9)	22.2 (13.1)	21.3 (12.6)	0.043
CRP (mg/L)	7.1 (6.9)	9.9 (9.8)	21.8 (18.9)	18.5 (21.5)	21.1 (21.1)	<0.001
HLA-B27-positive	50/54 (92.6%)	73/103 (70.9%)	72/82 (87.8%)	94/119 (79.0%)	44/55 (80.0%)	0.007
Presence of (bridging) syndesmophytes	16/39 (41.0%)	45/73 (61.6%)	38/70 (54.3%)	64/90 (71.1%)	18/34 (52.9%)	0.017
Start anti-TNF-α treatment	15/54 (27.8%)	41/101 (40.6%)	81/86 (94.2%)	85/120 (70.8%)	35/55 (63.6%)	<0.001
BASDAI (0–10)	1.7 (1.3)	4.0 (2.0)	6.0 (1.4)	5.9 (1.6)	7.3 (1.5)	<0.001
BASFI (0–10)	1.8 (1.8)	3.8 (2.3)	5.2 (2.1)	6.1 (1.8)	7.0 (1.9)	<0.001
Patient GDA (0–10)	1.8 (1.7)	4.3 (2.4)	6.9 (1.9)	6.3 (2.0)	7.3 (1.7)	<0.001
Physician GDA (0–10)	2.1 (2.0)	2.8 (1.9)	4.5 (2.2)	4.4 (2.3)	4.6 (2.3)	<0.001

Values represent mean (SD) or number (%).

*Number of patients per group at baseline.

anti-TNF-α treatment, antitumour necrosis factor α treatment; ASDAS-CRP, Ankylosing Spondylitis Disease Activity Score—C reactive protein; ASQoL, Ankylosing Spondylitis Quality of Life; BASDAI, Bath Ankylosing Spondylitis Disease Activity Index; BASFI, Bath Ankylosing Spondylitis Functional Index; BMI, body mass index; CRP, C reactive protein; GDA, Global Disease Activity Score; GLAS, Groningen Leeuwarden Ankylosing Spondylitis; HLA-B27, Human Leucocyte Antigen B27; HRQoL, health-related quality of life; OASIS, Outcome in Ankylosing Spondylitis International Study.

Overall, marked heterogeneity of disease impact on progression of HRQoL and disease activity was observed, which could be linked to both individual patients’ characteristics and treatment with TNF-α-inhibitors. In this respect, a complex picture emerged. Males with low CRP values, HLA-B27 positivity without (bridging) syndesmophytes, displayed favourable clinical manifestations, despite absence of TNF-α-inhibitors. At the same time, effectiveness of TNF-α-inhibitors on disease impact was observed, particularly in younger, male patients with shorter symptom duration. By contrast, older females with longer symptom duration were more prone to a disadvantageous clinical picture, showing only minor or no treatment responsiveness.

In both cohorts, and subsequently in each revealed trajectory group, losses to follow-up occurred. Disease impact appeared to be of no influence to dropout.

Discussion

The concept of HRQoL is an increasingly important health-related outcome, as it reveals where medical care can add value to patients’ lives. In axSpA, disease activity is the main modifiable factor contributing to HRQoL. While it is recognised that HRQoL and disease activity vary largely within as well as across patients, there is insufficient knowledge on how their temporal codependencies unfold. By means of bivariate GBTM, we were able to explore the development of HRQoL and disease activity simultaneously. In this study, five trajectories of HRQoL and disease activity joint evolution with different impact of disease have been identified: two trajectories with stable low/moderate impact (t1 and t2), two trajectories with relatively stable high/very high impact (t4 and t5) and one trajectory with major improvement in HRQoL and disease activity reaching low impact of disease (t3). These developmental trajectories provide several substantive insights: first, HRQoL and disease activity seem to move in parallel, as the ranking of both the HRQoL- and disease activity trajectories as well as their patterns of changes over time are to a great extent a reflection of each other; second, each group shows a relatively stable course of both HRQoL and disease activity, except for t4 and t3; third, the findings indicate that some patients may share similar perceived disease activity, but can still differ in the self-reported QoL (see t2 and t4 from time point 1 onwards, figure 1) and finally, each group shows some unique patient characteristics that may underlie the distinctiveness of their temporal courses.

The five revealed disease impact trajectories and their characterisations are largely consistent with the reported disease activity trajectories by Moltó et al.17 For example, the largest proportion of females, highest BASDAI-scores and BASFI-scores can be found in the trajectory very high impact (t5)/persistently very high disease activity. Conversely, younger patients in combination with prescription of TNF-α-inhibitors characterise the trajectory high impact with major improvement (t3)/changing disease activity. Although our study can be seen as both a confirmation of and a complement to the findings by Moltó et al, there are some important methodological differences. Moltó et al have analysed patients with early axSpA over 3 years, while in the present study, the majority of the patients had established (radiographic) axSpA and were followed up to 8 years, thereby expanding the knowledge on differences in disease-behaviour in those with longstanding disease. Given the similarities between our trajectories and the ones found by Moltó et al, it seems that the trajectories for patients with early disease continue in patients with established disease. Moreover, because of the longer follow-up period, we could also demonstrate that most disease activity trajectories remained, at the group level, stable over time. Only when treatment with TNF-α-inhibitors was initiated in patients with high disease activity, a (major) improvement in impact of disease was seen in most patients, after which the values stabilised again. Last, while Moltó evaluated the temporal course of disease activity only, we were able to capitalise on its coevolution with HRQoL, demonstrating, for the first time, the tight equivalence of their temporal unfolding.

For the GBTM analyses, data from two different cohorts, OASIS (2000–2008) and GLAS (2004–2012) were merged. Their timelines did not fully coincide, but both were followed biennially for up to 8 years. In the class-comparisons, the cohorts’ distribution emerged as significantly different among the latent trajectories. For example, t3 almost entirely existed of GLAS patients, while t1 is largely represented by OASIS patients. This quasicomplete separation hints at the presence of a cohort effect, that is, that the GBTM model clustered subjects based on their shared timeline. If so, identified clusters would likely represent an artefact explained by the impact of the groups bonded by distinct periods, and their common environments (exposures). This conclusion is plausible not because of an unknown period effect. The impact of the cohorts’ composition on class-extraction, we argue, is primarily attributable to confounding by inclusion criteria, which happened to be linked to the period. OASIS patients had an established diagnosis of axSpA at a time that biologicals were not yet prescribed (except for one patient included in a clinical trial). Conversely, nearly all GLAS patients were prescribed TNF-α-inhibitors because of active disease. Hence, the culprit underlying the differential cohort distribution can be pinpointed: their baseline clinical features, for example, disease duration, inflammatory components and associated treatment status. These factors are indisputably linked to the two outcomes. By means of example, the major improvement trajectory t3 was composed of young subjects with shorter disease duration, high level of inflammation markers and presence of biological. However, findings also raise the possibility of alternative explanations for data segmentation beyond the cohort effect. For instance, in one-third of the patients assigned to the very high impact trajectory t5 were OASIS members, who happen to display a similar clinical progression as their GLAS counterparts. Likewise, both cohorts were well represented in trajectory t2. More importantly, OASIS and GLAS did not differ in some features that emerged as discriminatory among the latent classes. Therefore, a cohort effect alone does not seem to fully account for the identified temporal patterns. An interplay of mechanisms underpinning the level and unfolding (eg, responsiveness to treatment) of the developmental paths is bound to offer a biologically more sound explanation for the findings.

Our study has several limitations. First, temporal precedence cannot be established with this design as both outcomes were measured simultaneously, that is, it cannot be concluded that disease activity changed before HRQoL or vice versa. Second, the extracted classes and their names are no true entities, but only features of the data, since each patient was assigned to one of the groups with a certain probability. We provided the names to help describe the trajectories in terms that would enhance a meaningful framework for clinicians and improve interpretability.

The main strengths of our study are the replication of findings of Moltó et al, the expansion of these findings by using an innovative methodology allowing evaluation of patterns of joint evolution of HRQoL and disease activity, and the opportunity to better clarify the role of anti-TNFα treatment on the impact of the disease. The main advantage of this bivariate model is that we could examine how these two variables simultaneously changed over time. The combination is expected to be more informative, besides being more attuned to clinicians’ demands of a more rounded clinical picture.

Our study has also several clinical implications. Our starting hypothesis was that the heterogenic clinical picture of axSpA might be further explored by assigning individuals to subgroups based on more than one parameter, and thereby improve diagnostic utility and clinical relevance of the identified data-driven classes. Patients in t1 showed an enduring low impact of disease with favourable HRQoL and low disease activity. This suggests that potentially, these patients may be followed-up with longer intervals between clinical visits. Similarly, patients in t3, which showed major improvement after initiation of anti-TNF-α treatment, with subsequently comparable low disease impact, can be followed-up with longer intervals. However, this needs to be demonstrated in separate studies. More burdensome are the patients with persistent high disease impact, despite anti-TNF-α treatment in some. Possibly, concomitant fibromyalgia plays a role in perceived disease activity and HRQoL. Previous studies showed that concomitant fibromyalgia may be present in approximately 20% of patients with axSpA, which is associated with poorer response to anti-TNF-α treatment.30 31 In our study, unfortunately no information was collected on a history of fibromyalgia. Nevertheless, it seems that this population would benefit from additional support for their impact of disease.

In conclusion, the five latent trajectory groups revealed the heterogeneous character of axSpA, the joint-evolution of HRQoL and disease activity, the different experiences of impact of axSpA among patients’ subgroups, and the characteristics of these subgroups. The evolution of HRQoL and disease activity patterns of change/stability are tightly linked, providing supportive evidence to a treat-to-target strategy envisaging improvement of patients’ HRQoL through management of disease activity.