Assessing the Course of Organ Dysfunction Using Joint Longitudinal and Time-to-Event Modeling in the Vasopressin and Septic Shock Trial.

Non-mortality septic shock outcomes (e.g., Sequential Organ Failure Assessment score) are important clinical endpoints in pivotal sepsis trials. However, comparisons of observed longitudinal non-mortality outcomes between study groups can be biased if death is unequal between study groups or is associated with an intervention (i.e., informative censoring). We compared the effects of vasopressin versus norepinephrine on the Sequential Organ Failure Assessment score in the Vasopressin and Septic Shock Trial to illustrate the use of joint modeling to help minimize potential bias from informative censoring.
DESIGN: Secondary analysis of the Vasopressin and Septic Shock Trial data.
SETTING: Twenty-seven ICUs in Canada, Australia, and United States.
SUBJECTS: Seven hundred sixty-three participants with septic shock who received blinded vasopressin (n = 389) or norepinephrine infusions (n = 374).
MEASUREMENTS AND MAIN RESULTS: Sequential Organ Failure Assessment scores were calculated daily until discharge, death, or day 28 after randomization. Mortality was numerically higher in the norepinephrine arm (28 d mortality of 39% vs 35%; p = 0.25), and there was a positive association between higher Sequential Organ Failure Assessment scores and patient mortality, characteristics that suggest a potential for bias from informative censoring of Sequential Organ Failure Assessment scores by death. The best-fitting joint longitudinal (i.e., linear mixed-effects model) and survival (i.e., Cox proportional hazards model for the time-to-death) model showed that norepinephrine was associated with a more rapid improvement in the total Sequential Organ Failure Assessment score through day 4, and then the daily Sequential Organ Failure Assessment scores converged and overlapped for the remainder of the study period.
CONCLUSIONS: Short-term reversal of organ dysfunction occurred more rapidly with norepinephrine compared with vasopressin, although differences between study arms did not persist after day 4. Joint models are an accessible methodology that could be used in critical care trials to assess the effects of interventions on the longitudinal progression of key outcomes (e.g., organ dysfunction, biomarkers, or quality of life) that may be informatively truncated by death or other censoring events.

The Sequential Organ Failure Assessment (SOFA) score (1) is frequently used to measure organ failure in sepsis randomized clinical trials (RCTs) (2–12) and is increasingly used as a primary trial outcome (2, 13). However, the use of an outcome such as the SOFA score presents difficult challenges for trial interpretation. For example, SOFA values may be worse for those who died before death and unmeasured (i.e., missing) after due to the competing risk of death. This type of informative censoring issue creates a missing data problem known as “missing not at random” (MNAR). As a result, comparisons of non-mortality outcomes between study groups can be biased and misleading, particularly if mortality is not equal (i.e., differential) between study groups.

Several approaches can be used to address the problem of informative censoring by competing risks in RCTs (14–16). Investigators can assign fixed values to those that die to create a composite endpoint such as event-free days (17), use an ordinal rank composite outcome such as the Rankin score (18), or a paired outcome ranking system such as the win ratio (19, 20). However, composite outcomes may lose details regarding the component outcomes (21, 22), which complicate their interpretation. So-called “joint models” are an alternative statistical framework that can be used to address the potential problems introduced by informative censoring (23–27), although their use in critical care research has thus far been limited (28–32). The aim of this article is to introduce the use of a joint longitudinal and time-to-event (survival) model approach to analyze critical care trials by evaluating the effects of vasopressin versus norepinephrine on the SOFA score through day 28 in the Vasopressin and Septic Shock Trial (VASST) (33). To support critical care researchers interested in applying this methodology, we have provided statistical code using the open-source R language (34) for statistical computing (Supplementary Digital Content, http://links.lww.com/CCX/A161).

METHODS

Study Design and Participants

This was an unplanned post hoc analysis of VASST (33), a multicenter, double-blind RCT that assigned patients who had septic shock and were receiving a minimum of 5 μg/min of norepinephrine to receive either low-dose vasopressin (0.01–0.03 U/min) or norepinephrine (5–15 μg/min) plus open-label vasopressors. Written informed consent was obtained from all participants. The analysis was approved by the Institutional Review Board at the University of British Columbia.

SOFA Score

The total SOFA score was calculated from the six organ-specific subscores using data recorded in the VASST case report form. Participants with at least one calendar day in which all SOFA subscores were recorded were included in analysis. Intermittent missingness of SOFA subscores before discharge or death was considered to be missing at random (MAR). Missing values after discharge or death were considered MNAR and thus nonignorable.

Rationale for Using a Joint Model

Comparisons of non-mortality longitudinal endpoints (e.g., SOFA) are susceptible to biases from informative censoring or truncation due to death or discharge, or other competing risks that occur during the observation period (14, 35–37). A traditional mixed-effects model is robust to outcomes that are MAR, whereas a joint model can handle data that is MNAR (e.g., from informative censoring) (38, 39). Joint modeling entails the simultaneous estimation of two separate regression models with a shared random effect (see References [23, 24] for a technical tutorial). The most common joint model combination includes a longitudinal (also called “repeated measures” or “mixed-effects models”) and time-to-event (survival) model. The underlying models retain their familiar interpretations. However, the “simultaneous” or “joint” modeling of the longitudinal and survival processes allows each model to inform (i.e., adjust) the other. The sharing of information between the two models can help mitigate potential biases caused by missing outcome data in the longitudinal model (i.e., SOFA score in this study) because the post-randomization event causing missing data is explicitly incorporated in the joint model. This approach to modeling the longitudinal SOFA outcome entails the estimation of an unobserved, or latent, SOFA trajectory that provides an estimate of how the SOFA trajectory may have looked if there was no informative censoring (e.g., death), and all observations were observed.

Statistical Analysis

The focus of our analyses was the comparison of the effects of vasopressin versus norepinephrine on the total SOFA score. To assess the impact of mortality on the treatment effect of vasopressin versus norepinephrine on SOFA score, we fitted a traditional linear mixed-effects model of the SOFA score over time with the same model formulation of the best-fitting joint model for comparison purposes, and several joint longitudinal and survival models using the formulation of Henderson et al (40). SOFA scores for participants were available through their day of hospital discharge, death, or day 28, whichever occurred first. The time-to-event survival data were modeled using a cause-specific Cox proportional hazards model for death. Hence, both discharge and day 28 were considered as censoring events in the death-specific Cox model. Shared random effects were used to capture the association between the longitudinal and time-to-event submodels.

We explored multiple model specifications for the joint model by varying the functional form of time in both the fixed effects in the longitudinal component and the shared random-effects specifications. For instance, we tested linear and quadratic effects, and spline terms with an increasing number of knots (up to 7 and 2 for the fixed effect and random effect of time, respectively). We selected the final model according to the Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC). ses for the coefficients of the joint model were estimated via 1,000 bootstrap replications with resampling at the participant level. We tested for the modification of treatment by study time by including an interaction term in the longitudinal submodel and by testing its significance using the Wald test. Finally, we performed several sensitivity analyses to assess the robustness of our findings (Supplementary Digital Content, http://links.lww.com/CCX/A161), including fitting a joint model with a cause-specific model for discharge and a model for either death or discharge.

All models were estimated using the joineRML package (version 0.4.2 [41]) in R (34). Joint models for longitudinal and survival data can also be implemented in other statistical software such as Stata (StataCorp LLC, College Station, TX) (42) and SAS (SAS Institute, Inc., Cary, NC) (43).

RESULTS

In VASST, 779 participants were randomized (397 received vasopressin and 382 norepinephrine). We included 763 participants (97.9%, n = 389 vasopressin, n = 374 norepinephrine) who had sufficient data to calculate at least one total SOFA score in our analytic sample; eight randomized participants per arm were excluded who did not have a measured SOFA score. Our analytic sample included 6,934 SOFA score measurements (n = 3,476 vasopressin, n = 3,458 norepinephrine; Fig. ), with a median of seven SOFA measurements per participant (interquartile interval [IQI], 4–12).

Proportion (A) and number (B) of individuals within each Sequential Organ Failure Assessment (SOFA) score subgroup over time in the Vasopressin and Septic Shock Trial. Missing values are included and labeled as not applicable (NA).

Mean daily SOFA scores throughout the study period are shown in Figure by survival status and in Figure by study arm. As the study period progressed, more SOFA scores were missing or censored, resulting in fewer observations per participant toward the end of the study period. Patients who died had higher SOFA scores compared with those who survived (Figs. 2, A and B; median of 13, IQI, 10–15 vs a median of 11, IQI, 8–13; Mann-Whitney U test p < 0.01), and patients randomized to receive vasopressin had slower early declines in SOFA than patients receiving norepinephrine alone (Fig. 3A), appearing to be driven by differences in the cardiovascular SOFA subscore (Fig. 3B). In unadjusted analysis, there was a numerically higher 28-day mortality observed in the norepinephrine group (35.0% vs 39.3%; χ2 p = 0.25) (Fig. ).

Figure 2.

Progression of the average total Sequential Organ Failure Assessment (SOFA) score (A) and the SOFA organ subscores (B), overtime among those who survived and died in the Vasopressin and Septic Shock Trial. Survivors are marked by the solid line, and nonsurvivors are marked by the dashed line.

Figure 3.

Progression of the average total Sequential Organ Failure Assessment (SOFA) score (A) and the SOFA organ subscales (B), overtime by treatment arm in the Vasopressin and Septic Shock Trial. The vasopressin group is marked by the solid line, and the dashed line denotes the norepinephrine group.

Joint Longitudinal and Survival Analysis

We examined 36 potential joint model specifications (Supplementary Digital Content, http://links.lww.com/CCX/A161)and selected the best-fitting model using AIC and BIC (Fig. ). The model with the best fit had the fixed effect of time modeled using a natural spline with 7 degrees of freedom, a random intercept, and a random effect of time modeled using a natural spline with 2 degrees of freedom.

The best-fitting model did not demonstrate an association between vasopressin versus norepinephrine treatment and the survival component of the joint model, consistent with the original findings in the primary trial report that vasopressin did not decrease mortality (hazard ratio of norepinephrine vs vasopressin: 1.15; 95% CI, 0.85–1.46). However, there was a significant association between the longitudinal SOFA score and the risk of death (association parameter = 0.25; z test p < 0.01). This association indicates that changes in the SOFA score over time were associated with 28-day mortality. The estimated coefficients from the joint model are reported in the Supplementary Digital Content (http://links.lww.com/CCX/A161).

The joint model showed evidence of an interaction between time and treatment, suggesting that the longitudinal evolution of SOFA was different between the vasopressin and norepinephrine groups over time (Fig. 5A) (Wald test χ2 = 42.0; p < 0.01). Specifically, the norepinephrine arm continued to show a larger decrease in the SOFA score than in the vasopressin arm at day 2 (SOFAVasopressin – SOFANorepinephrine = 1.17; 95% CI, 0.46–1.87) and day 3 (1.20; 0.42–1.99). However, starting on day 4, the magnitude of the difference in SOFA scores (0.70; –0.03 to 1.42) began to diminish and then overlapped between arms through the end of the follow-up period on day 28 (Fig. 5). Similar results were found using a traditional linear mixed-effects analysis (Fig. 5; e.g., day 2: SOFAVasopressin – SOFANorepinephrine = 1.14; 95% CI, 0.70–1.59; and day 3: 1.16; 95% CI, 0.70–1.63.)

Figure 5.

Model of the total Sequential Organ Failure Assessment (SOFA) score over time using a joint model (JM) and mixed-effects (MEs) model in the Vasopressin and Septic Shock Trial. The longitudinal trajectories for both models are presented in (A), and the difference over time between the vasopressin and norepinephrine treatment arms is presented in (B). The vasopressin group is marked by the solid line, and the norepinephrine group is marked by the dashed line. ses are estimated using 1,000 bootstrap replications. The joint model used a natural spline with 7 degrees of freedom (df) for the fixed effect, a random intercept, and a random effect of time modeled using a natural spline with 2 df. The joint model took 1.7 min to fit without accounting for the 1,000 bootstrap replications.

Finally, we fit the cause-specific model for the competing risk of discharge using the same formulation for the fixed effects and random effects of the best model selected by the AIC and BIC. The estimated longitudinal SOFA score trajectory was similar to the estimated trajectory from the main analysis model, and the treatment was not associated with the risk of discharge (log hazard ratio, –0.08; 95% CI, –0.40 to 0.24). Our results were also stable when using an any-event (discharge or death) model. A detailed technical summary, including statistical code for replication, additional data summaries, along with a description of joint models in general, and a description specific to our final model, is provided in the Supplementary Digital Content (http://links.lww.com/CCX/A161).

DISCUSSION

In this post hoc secondary analysis of the VASST trial, we used joint linear mixed-effects and Cox proportional hazards models to compare the trajectory of SOFA scores over time between the vasopressin and norepinephrine study arms. Despite a small differential occurrence and timing of death between treatment arms (Fig. 4), the findings of joint modeling and conventional modeling with mixed-effects regression were similar (Fig. 5). Both approaches demonstrated that SOFA scores improved more quickly among patients receiving norepinephrine when compared with vasopressin; however, this did not correlate with differences in mortality between arms. Although our use of joint models did not substantively change the interpretation of the VASST trial, we present our findings to demonstrate the potential utility of joint models in future trials in which a larger difference in the competing risk of death (or other censoring events) between trial arms may produce a stronger bias in non-mortality endpoint comparisons.

Figure 4.

Kaplan-Meier survival curves comparing study participants in the vasopressin group and norepinephrine group in the Vasopressin and Septic Shock Trial.

There is no ideal strategy for dealing with the competing risk of death in critical care trials when a non-mortality outcome is of interest. Unless an intervention truly has no impact on mortality, researchers must consider the potential for bias in non-mortality treatment effect estimates. For example, in a simulation study of methods to compare length of stay, we observed that even small mortality differences (i.e., 2.5% and 5% on the absolute scale) can lead to biased comparisons when using several common statistical models (15). Presently, many trialists deal with the occurrence of death by creating composite endpoints (17, 44). Alternative statistical modeling proposals put forth in the critical care literature include the gamma mixture model (45), the Fine and Gray model (46), principal stratification (47, 48), and the focus of this article, joint longitudinal and time-to-event (or survival) models (28).

We believe that joint models offer several inferential benefits to researchers. Foremost, the underlying longitudinal and survival models retain their well-understood interpretations. Although we used a linear mixed-effects and Cox proportional hazards model in our analysis, the joint modeling framework could be extended to any longitudinal clinical outcome (e.g., blood pressure, fluid balance, urine output, arterial pH, arterial lactate concentration, or daily presence of delirium [28]) and more complex survival models. Different distributions and assumptions would simply use alternate longitudinal (e.g., binary or count data) and/or time-to-event (e.g., competing risks) specifications. Second, adjustment for prognostic variables and center effects can be straightforwardly implemented. Another benefit of the joint model is the adjustment of the mortality analysis to account for the values of the longitudinal outcome (e.g., severity of organ dysfunction) over time (i.e., time-varying covariate adjustment [23, 24]). The use of a joint model may offer particular interpretive benefits over composite endpoints or event-free day measures (e.g., organ failure-free days) (17, 44). For example, although composite endpoints may capture the “net effect” of an intervention, they reduce detailed longitudinal data into a single value, losing benefits that denser and more detailed longitudinal data provide about the trajectory and modification of longitudinal health states due to an intervention (45).

There are limitations and challenges to using the proposed joint modeling approach. First, joint models can be computationally intensive. However, statistical computing in this area is advancing rapidly, making these models increasingly accessible in common software packages. Second, a simple numerical summary of a difference over time may not be straightforward if trajectories are nonlinear, as observed in our study (Fig. 5). Thus, we relied on the assessment of daily differences and visual depictions to summarize our results. Although these strategies can produce clinically informative knowledge, they cannot produce simple summaries that may be desired in some scenarios, such as in a regulatory submission. Third, for illustration, we estimated cause-specific survival models that assumed no competing risks. Thus, survival (or discharge) probabilities estimated from our model would not be informative. Although the use of a cause-specific model does not invalidate our results given the focus on SOFA trajectories overtime, there are several potential questions in critical care where a competing risks specification may be preferable.

Our analysis identified important questions that need to be addressed to support the use of SOFA scores as an outcome. Foremost, there is a need for consensus on a clinically relevant change or difference in SOFA scores between intervention arms. Second, it is unclear how an acute change or difference in SOFA scores between intervention arms translates into long-term non-mortality outcomes (e.g., quality of life). Although we did not observe a short-term mortality benefit associated with SOFA declines in the current study, researchers have shown an association between lower short-term organ dysfunction and improved long-term survival, suggesting that sepsis therapies that reduce short-term organ dysfunction may lead to better long-term outcomes (49, 50). The optimal integration of the SOFA score into critical care trials is an important topic for future research, especially given its use as a primary outcome in recent trials (2, 13).

An additional limitation of our study is that we did not examine the SOFA subscores in the joint model regression framework because the small number of categories (0–4) does not lend itself to straightforward or informative modeling strategies. However, Figure 3B suggests that the cardiovascular SOFA score recovery was slower in the vasopressin arm compared with the norepinephrine arm, and this might be what we observed in the overall SOFA score analysis. The calculated vasopressor dose in VASST included both open-label and blinded study drug. Differences in cardiovascular SOFA are not likely due to differences in drug potency because vasopressors were titrated to a mean arterial pressure target in both groups. Differences in early cardiovascular SOFA are also not likely caused by differences in early rates of death (survival curves show overlap of vasopressin and norepinephrine groups until day 10). As such, the explanation for this organ-specific subscore difference is unclear.

To conclude, this reanalysis of the VASST trial provided a case study of how joint longitudinal and survival models could be used to augment assessments of non-mortality outcomes at risk for informative censoring in critical care research. Using the joint modeling framework, we modeled complex, nonlinear relationships of the SOFA score over the course of the VASST trial. Our best-fitting model suggested a slight benefit in improvement of organ dysfunction (as measured by the SOFA score) in the norepinephrine arm compared with vasopressin in the first few days of the study, but this change was not associated with differences in mortality. Although joint modeling did not produce substantially different estimates than traditional modeling methods in this study, we demonstrate the application of joint models in critical care trials. It will be important that researchers who use joint models report their model structure, results, and interpretation clearly and accurately in order to subsequently use the results in future validation studies and evidence synthesis (51).

ACKNOWLEDGMENTS

We thank the patients and families who participated in the Vasopressin and Septic Shock Trial study which has increased our understanding of the treatment of septic shock. We also thank the many dedicated clinicians (doctors, nurses, therapists and others) who cared for these critically ill patients and comforted their families.