Defining phenotypes and treatment effect heterogeneity to inform acute respiratory distress syndrome and sepsis trials: secondary analyses of three RCTs

Conceptual approach for heterogeneity of treatment effect

Our aim was to assess whether or not an individual patient’s baseline risk of death modifies the treatment effect of an intervention (HTE). The Acute Physiology And Chronic Health Evaluation II (APACHE II) model has been proposed as a potential model for HTE evaluation. 24,33,34 We assessed HTE using the APACHE II score34 as the primary measure of baseline risk. In addition, we assessed three secondary measures based on the APACHE II model: (1) the APACHE II physiology score (APS-APII), (2) the APACHE II-calculated risk of death (R_calc. ), as originally proposed by Knaus and colleagues,34 and (3) a modified APACHE II risk of death model recalibrated (R_recal. ) using data from the VANISH trial15 and the LeoPARDS trial. 16 The rationale for using the APS-APII was that the total APACHE II score determines a non-modifiable risk of death based on age and severe comorbidity, but the physiological derangement most likely mediates the relationship between treatment effect and outcome. 35 We also investigated whether or not any HTE could be driven by adverse events: if low-risk patients have similar exposure to treatment-related harms as high-risk patients, but do not have the same exposure to benefits, this would result in a net harm signal. 24 Furthermore, irrespective of whether the treatment effects of interventions varied or remained constant over the range of baseline risk, HTE may manifest because of differences in treatment-related adverse events over the range of baseline risk.

The VANISH trial

The VANISH trial15 is a 2 × 2 factorial, double-blind RCT in adult patients with sepsis who required vasopressors and was carried out in 18 general adult ICUs in the UK. In the VANISH trial, patients were randomly allocated to vasopressin and hydrocortisone (n = 101), vasopressin and placebo (n = 104), noradrenaline and hydrocortisone (n = 101) or noradrenaline and placebo (n = 103). Patients received the second study drug (i.e. hydrocortisone/placebo) only if the maximum infusion of the first study drug (i.e. vasopressin/noradrenaline) had been reached. The 28-day mortality was 63 of 204 (30.9%) patients in the vasopressin group and 56 of 204 (27.5%) patients in the noradrenaline group [a difference of 3.4%, 95% confidence interval (CI) –5.4% to 12.3%].

The LeoPARDS trial

The LeoPARDS trial16 is a two-arm, parallel-group, double-blind, placebo-controlled RCT in adult patients with sepsis who required vasopressors, carried out in 34 ICUs in the UK. In the LeoPARDS trial, patients were randomised to receive either levosimendan (n = 258) or placebo (n = 257) over 24 hours, in addition to standard care. The 28-day mortality was 89 of 258 (34.5%) patients in the levosimendan group and 79 of 256 (30.9%) patients in the placebo group (a difference of 3.6%, 95% CI −4.5% to 11.7%).

The HARP-2 trial

The HARP-2 trial17 is a two-arm, parallel-group, double-blind, placebo-controlled RCT in adult patients (within 48 hours after the onset of ARDS), carried out in 40 ICUs in the UK and Ireland. In the HARP-2 trial, patients were randomised to receive either once-daily simvastatin or identical placebo tablets enterally for up to 28 days. The 28-day mortality was 57 of 259 (22.0%) patients in the simvastatin group and 75 of 280 (26.8%) patients in the placebo group [risk ratio (RR) 0.8, 95% CI 0.6 to 1.1].

Groups for comparison

Treatment effects were assessed primarily on an intention-to-treat basis, except for hydrocortisone compared with placebo (i.e. the second comparison in the VANISH trial15). Patients were eligible to receive hydrocortisone/placebo only if they had reached the maximum infusion of the first study drug, which occurred for around three-quarters of patients. This eligibility criterion was applied post randomisation, but before the administration of the second (blinded) study drug. As there was no interaction between the study drugs, and given the limited power of the analysis, only patients eligible to receive the second drug were included in this comparison (hydrocortisone, n = 148; placebo, n = 148). A sensitivity analysis on the per-protocol population was conducted for the other drug comparisons. Each trial was analysed separately.

Outcomes

The primary outcome is mortality at 28 days after randomisation. The secondary outcome is hospital mortality during the initial hospital stay (i.e. ignoring readmissions). The rationale for this is that patients in these trials who left hospital alive were either well enough to be discharged or still sick but transferred elsewhere (e.g. social care). Therefore, hospital mortality is not a true binary outcome, as those alive at discharge are not a consistent group; hence landmark mortality at 28 days was preferred.

Measures of baseline risk

The primary analysis examined HTE for 28-day mortality, with the APACHE II34 as the measure of baseline risk, comparing treatment effect in patients with an APACHE II score above (high) or below (low) the overall median score of 25 points. This score has already been suggested a measure over which HTE could be evaluated. 24,33 The APACHE II score is the sum of the points from three elements: (1) acute physiology, (2) age and (3) chronic health. The calculation of APACHE II score is given in Table 1.

As secondary analyses we examined three other baseline risk measures. The first is the acute physiology element of the APACHE II score, which we denote APS-APII. The rationale for using the APS-APII was that the total APACHE II score includes non-modifiable risk of death attributable to age and comorbidity, but the physiological components are more likely mediators of the effect of treatment on outcome. 35 The second additional baseline risk measure we considered was the risk of death in hospital (i.e. R_calc. ), calculated based on APACHE II score using a formula originally proposed by Knaus and colleagues34 as follows:

Post-emergency surgery is a binary indicator, and diagnostic category weights relate to the principal reason for admission for a patient. The third baseline risk measure was R_recal. , given by recalibrating R (risk) to the study population to see whether or not an improved prediction yielded a different estimate of HTE. Methods for developing R_recal. are given in Recalibrating APACHE II.

Recalibrating APACHE II

The following logistic regression models were estimated, with each subsequent model recalibrating with finer detail:

where AP is acute physiology and AP_{1, . . . , 12} is temperature, mean arterial pressure, heart rate, respiratory rate, oxygenation, arterial pH, serum sodium concentration, serum potassium concentration, serum creatinine concentration, haematocrit, white blood cell count and Glasgow Coma Scale score, respectively, and CH is chronic health. These elements were entered into the model on the points scale described in Table 1, as were age and chronic health. Mortality at 28 days was used as the outcome, as this was the outcome of interest, although predictive performance for hospital mortality was also assessed. Three additional models were estimated by adding the number of organ dysfunctions at baseline (respiratory, renal, hepatic, coagulation and cardiovascular), based on a SOFA score of ≥ 2 points. 7 The number of organ dysfunctions was treated as a continuous variable. The coefficients for the diagnostic categories were kept in proportion to the existing weights, rather than re-estimating the weights for each category because of sparse data.

The discriminatory performance of the models for the whole cohort was compared using the area under receiver operating characteristic (AUROC) curve and the discrimination slope (DS), which is the mean difference in prediction comparing those with the event and those without. 45 We did not use Cox recalibration [i.e. the prediction resulting from a logistic regression of the outcome of interest against the logit (R)], as this would not change the discrimination of the model. If the patients were ordered with respect to their score after Cox recalibration then their rank would remain the same and therefore the HTE pattern would also be the same.

The AUROC curve was calculated for both hospital and 28-day mortality. As some of the data on which performance is assessed were also used to build the model, bootstrapping was used to correct the AUROC curve for overoptimism. 46 For this method, a bootstrap sample is taken and the model is estimated on the sample to obtain new coefficients. In addition, the AUROC curve (for example) is calculated (AUROC_boot. ).

Next, the same model and coefficients are applied to the original data set and the AUROC curve calculated (AUROC_orig. ). Then, AUROC_boot.  – AUROC_orig. gives an estimate of the optimism and this is repeated for many bootstrap samples and the average optimism taken. The averaged optimism is subtracted from the optimistic AUROC to give a corrected AUROC. This process was modified because the model was estimated on only a sample of the data. A bootstrap sample of the whole data set was taken, stratifying on treatment (control vs. active, taking any active treatment for the VANISH trial15) to ensure that the proportion of placebos is the same in the bootstrap sample. AUROC_boot. is calculated by estimating the model on the placebo groups and applying it to the whole bootstrap sample and AUROC_orig. is calculated by applying the same model to the original data set. The same approach was applied to the DS. The model with the best corrected discriminatory performance was used as an additional measure of baseline risk.

Models were estimated using the control groups from the VANISH15 (noradrenaline + placebo, n = 103) and LeoPARDS16 trials (placebo, n = 257) to avoid using post-randomisation outcomes to calculate baseline risk. However, in a simulation study, Burke and colleagues47 found that using the whole cohort slightly reduced bias, overfitting and risk of a false-positive finding for HTE, and so in addition the recalibration was performed using the whole cohort.

Descriptive analysis

Distributions of the baseline risk measures in the trial populations were described with histograms, including by treatment group, to check whether or not the distribution was balanced. APACHE II score was grouped in increments of 5 points, with those scoring ≥ 35 points in one category (this is the same categories used by Knaus and colleagues34) and risk of death was grouped into 10% increments. The relationship between risk measures and mortality in the trial cohorts was described using bar charts showing the proportion of patients who died in each category. The discriminatory performance was assessed using the AUROC curve. We estimated the extreme quartile odds ratio (EQuOR) (i.e. the ratio of the odds of death in the highest vs. lowest quartile for risk) as an estimate of how the risk of death varies between patients in the same trial. 48

Regression modelling

We also considered heterogeneity over the range of APACHE II as a continuous variable. A logistic regression model was constructed with 28-day mortality as the outcome, and treatment, APACHE II score and an interaction between the two as covariates to test whether or not treatment effect varies over baseline risk. Logistic regression was chosen over Cox regression as the outcome is short term (i.e. 28 days) and there was no censoring. The model is given by:

where p is the predicted probability of death before 28 days, risk is the mean-centred APACHE II and trt is the binary treatment indicator. The effect of a unit increase in baseline risk on mortality in the placebo group is given by β_R and the treatment effect for someone with mean baseline risk is given by β_T. HTE is described by the interaction term (i.e. the additional treatment effect for each unit increase in baseline risk) denoted by γ [i.e. all effects are log-odds ratios (ORs)]. Relative HTE was quantified by the interaction between APACHE II score and treatment, expressed as a ratio of ORs. Additive HTE was illustrated by plotting the estimated absolute difference in mortality between treatment groups across the range of APACHE II.

In the first instance, APACHE II was treated linearly, with residual plots used to determine if any transformations or non-linear terms were necessary. This would occur if, for example, the high-risk patients are too sick to benefit from the treatment, resulting in a n-shaped pattern of treatment effect. The non-linearity of the effects of baseline risk and of the interaction was investigated by grouping risk measures into quartiles and comparing nested models with linear and categorical associations using likelihood ratio tests.

Sensitivity analyses

Two sensitivity analyses for the main baseline risk measure (APACHE II score) were performed. First, we used hospital mortality as the outcome instead of mortality at 28 days, as the APACHE II score was originally devised as a prediction tool for hospital mortality. Second, we investigated the potential impact of missing data on the results. There were 47 patients in the VANISH trial15 for whom at least one element of the acute physiology score was missing (and 61 patients in the LeoPARDS trial16). In the main analysis, normal scores were assumed for these elements, as for the main trial. Total APACHE II scores were missing for 66 patients in the HARP-2 trial,17 and these patients were omitted from the main analysis, but are included in the forest plot. Missingness occurred pre randomisation and hence is independent of treatment effect; however, it may affect the precision of the results. In the sensitivity analysis we assumed that patients with missing data were (1) equally likely to be in the high-risk group as those with complete data, (2) 10% more likely or (3) 10% less likely. The APACHE II category was imputed 20 times under these assumptions, and the difference in RD and ratio of RR was computed as for the main analysis, combining results across imputations using Rubin’s rules. 49

Biomarker measurements

We measured three groups of markers to help delineate specific biological effects and illness characteristics. First, a limited cytokine profile was carried out to assess the balance between pro- and anti-inflammatory states using interleukins [i.e. interleukin 1 beta (IL-1β), IL-6, IL-8, IL-10, IL-17 and IL-18], soluble tumour necrosis factor receptor 1 (sTNFR1) and C–C motif chemokine ligand 2 (CCL2). The state of neutrophil and endothelial injury was assessed using myeloperoxidase (MPO), soluble intracellular adhesion molecule and angiotensin II (ANG II). For organ dysfunction, in addition to the SOFA variables,50 which were collected as part of trial data, we measured troponin and N-terminal pro-B-type natriuretic peptide (NT-proBNP) for cardiac dysfunction. These measurements used enzyme-linked immunosorbent assay-based methods. We had laboratory-specific standard operating procedures for these measurements prior to starting any measurements.

Exploratory analysis

Histograms and pairwise correlations were used to assess distributions, outliers and skewness (highly likely for the cytokine data). For assay data, the number and percentage of values below or above the limits of detection were recorded. Normal distributions were used for the continuous indicators, applying natural log transformations as necessary. Observations above or below the limits of detection were included in the analysis but treated as censored and all variables were standardised to have a mean of 0 and standard deviation (SD) of 1, with parameters taken from the data within the limits of detection. The number and proportion of missing observations were described for all biomarkers, clinical variables and demographic characteristics. If a patient has any missing individual indicators then LCA still allows the rest of the complete data to be included, implicitly assuming that the data are missing at random (i.e. the probability of missingness depends on only the observed data and not any missing data). This is reasonable for the biomarker data, as missing individual indicators are likely to be due to a technical issue.

Latent class modelling

Analysis was carried out separately for the LeoPARDS trial16 and VANISH trial15 cohorts. All measured biomarkers [i.e. PaO₂/FiO₂ ratio, creatinine, platelets, bilirubin, IL-1β, IL-6, IL-8, IL-10, IL-17, IL-18, MPO, soluble intercellular adhesion molecule (sICAM), ANG II, troponin, NT-proBNP, sTNFR1, lactate and CCL2] were included as indicator variables characterising the latent classes. Other baseline clinical and demographic variables {i.e. age, ethnicity, body mass index (BMI), comorbidities [any of New York Heart Association class IV (NYHA IV), severe chronic obstructive pulmonary disease (COPD), chronic renal failure, cirrhosis, immunodeficiency], site of infection (i.e. lung, abdomen, urine, other), SOFA score, APS-APII and post-surgical admission} that may be predictive of subphenotype were included in the model as class predictors. We also included APS-APII as a covariate in the submodel for each indicator based on a priori expectation of associations within classes. Only pre-randomisation data were used to develop the latent class model. All biomarkers were log-transformed and standardised because of skewness. Observations outside the limits of detection were included but treated as censored.

Latent class analysis models were fitted in three stages. First, conditional independence was assumed (i.e. all covariances constrained to zero) and no covariates predicting class membership were included. Second, prespecified clinical and demographic variables measured at baseline were included as covariates predicting class membership. Third, variance assumptions concerning indicators were relaxed to allow (1) non-constant residual variance across classes, (2) non-zero covariances and (3) both of these. It was not possible within the software used to model covariances between censored variables. For each stage we first fitted a one-class model, and then increased the number of classes by 1 until convergence could not be achieved. A number of strategies were used to achieve convergence, namely (1) for a k-class model, using starting values from a k – 1 class model; (2) using alternative integration methods; (3) reducing the number of censored indicators by treating values outside the limits as having values equal to the limit, for indicators with fewer than five such values; and (4) reducing the number of class predictors, selecting covariates that improved model fit based on likelihood ratio tests. Models were fitted using the gsem package in Stata® 15 (StataCorp LP, College Station, TX, USA).

The class means were estimated for each LCA model and differences across classes compared to determine which indicators showed the most separation across classes. For each model and each participant, the probability of an individual being in each class is predicted, with the probabilities for a participant summing to 1 across the classes. Each participant can then be assigned to the class for which they have the highest class probability. The Bayesian information criterion (BIC) was the primary measure of model selection,1 with smaller values indicating better fit. We also considered the Akaike information criterion (AIC), log-likelihood, entropy (i.e. a measure of class separation between 0 and 1), class sizes (with very small classes being indicative of overfitting) and the mean probability of class assignment, averaged over participants in the class. 51 We also assessed the class means and sized to see if the substantive interpretation of the classes differed across models. Additionally, plots of the change in fit statistics with the number of classes were used to determine where additional classes gave limited improvement in fit. 52 If models of different complexity gave a similar fit, then the simplest model was favoured.

Non-technical description of latent class analysis methods

Latent class analysis of the baseline variables aimed to replicate the previous publications using data from published ARDS trials. 23,42 We used data variables from subjects in all trial arms, without the influence of arm. The baseline variables consisted of clinical data, cytokine, and epithelial and endothelial injury marker profiles. For ARDS, the resulting subphenotypes were compared with the two subphenotypes derived independently in three trials previously. 23,42

The inclusion of variables, and any adaptation to their form, will depend on their robustness for their multivariate purpose, which was assessed by screening the univariate and bivariate data distributions for influential outliers, marked skewness and multicollinearity, for categorical variables with extreme prevalence and for variables contributing to the accumulation of missing data. This led to establishing the principal data set for the LCA of each trial, where the variables are further standardised to the z-scale to have mean of zero and unit variance, accounting for their differing units of measurement.

The latent class modelling stage involved the estimation of linear combinations of the standardised variables to identify a number of underlying classes. The number of classes will be determined formally by using the BIC and other model selection criteria, and by assessing the clinical interpretability of the classes as subphenotypes. With high probabilities of class membership, participants were assigned to their most likely phenotype. Regression methods with likelihood ratio tests were used to assess the association of classes with clinical outcomes, with randomisation kept intact and extended to compare response among randomised treatments. Given the factorial nature of the VANISH trial,15 this will involve a sequence of interactions tests respecting the design.

Description of subphenotypes

Once the most suitable model was selected, the estimated class means of each standardised indicator and their relative importance in class separation were shown by plotting the means, ordered by the magnitude of the largest difference between classes. Trial participants were assigned to the class for which they had the highest posterior probability of class membership for subsequent analysis. The median and interquartile range (IQR) of the observed biomarker values by class were tabulated, along with baseline clinical characteristics.

Clinical outcomes

For this study the primary clinical outcome was survival at 3 months in the LeoPARDS trial16 cohort, as this is the time point at which treatment differences stabilise. 53 Mean total SOFA score over 28 days (or ICU stay, whichever is shorter), which was the primary outcome in the LeoPARDS trial,16 and survival to 28 days were examined as secondary outcomes. For the VANISH trial15 we examined survival to 28 days (as survival to 3 months was not available), survival free of renal failure to 28 days among patients not in renal failure at baseline, and days alive and free of renal failure up to 28 days for all other patients (i.e. those who died or experienced some renal failure by day 28).

All outcomes were first compared between classes, irrespective of treatment, then treatment differences were compared between classes. For binary outcomes we presented the proportion of patients having the event in each class and performed a chi-squared test for the difference across classes. Treatment effects were expressed as a RD and the difference in treatment effects across classes as the a difference in RD. Ninety-five per cent CIs for the RD and difference in RD were calculated using linear regression with robust SEs. 54 For mean total SOFA score we presented the mean and SD, with differences between classes or treatment arms expressed as a difference in means. As mean total SOFA score is skewed, 95% CIs were calculated with bootstrapping, as was done in the main trial analysis. The median and IQR was presented for days alive and free of renal failure, again with bootstrap CIs. For continuous variables permutation tests were used to calculate p-values for the treatment–class interaction. Treatment effects by class were displayed using forest plots.

For the LeoPARDS trial,16 the first trial we analysed for identifying subphenotypes, we constructed a model to predict latent class, using a reduced set of indicators. A series of multinomial logit models were estimated, with latent class as the outcome and an increasing number of biomarkers as predictors, added in the order of greatest separation between classes. The probability of being in each class was predicted for each patient and patients were assigned to the class with the highest probability (similarly to the LCA). The class-specific sensitivity, specificity and c-statistics for each model were calculated by comparing the ‘gold-standard’ class of the latent class model with the ‘test’ class of the multinomial model. The final number of markers was chosen as the model for which the addition of further variables would bring negligible increases in accuracy measures.

Sensitivity analysis

In the main analysis we drew a distinction between class-defining and class-predicting variables. As a sensitivity analysis, we compared the class groupings when including all variables as indicators in the latent class model, following earlier work by Calfee and colleagues. 23,31

Descriptive analysis

In the VANISH,15 LeoPARDS16 and HARP-217 trials, 28-day mortality was not significantly different between the intervention and control arms (see Table 1). The illness severity (using the total APACHE II score) was lower in the HARP-2 trial17 than in the VANISH15 and LeoPARDS16 trials (see Table 1). The EQuOR highlighted significant heterogeneity of risk of death in all three RCTs for all three risk measures. Trial-level summary characteristics are shown in Table 2.

Figures 1–4 show the distribution of APACHE II and R by treatment arm in the VANISH trial15 (vasopressin vs. noradrenaline and hydrocortisone vs. placebo, shown separately). Figures 5 and 6 show the distribution of APACHE II and R, respectively, by treatment arm in the LeoPARDS trial16 cohort. The distribution of APACHE II in the HARP-2 trial17 cohort is shown in Figure 7.

FIGURE 1.

Distribution of APACHE II score by study drug 1: the VANISH trial15 cohort. (a) Noradrenaline; and (b) vasopressin.

FIGURE 2.

Distribution of APACHE II score by study drug 2: the VANISH trial15 cohort. (a) Placebo; and (b) hydrocortisone.

FIGURE 3.

Distribution of R_calc. by study drug 1: the VANISH trial15 cohort. (a) Noradrenaline; and (b) vasopressin.

FIGURE 4.

Distribution of R_calc. by study drug 2: the VANISH trial15 cohort. (a) Placebo; and (b) hydrocortisone.

FIGURE 5.

Distribution of APACHE II score by treatment: the LeoPARDS trial16 cohort. (a) Placebo; and (b) levosimendan.

FIGURE 6.

Distribution of R_calc. by treatment: the LeoPARDS trial16 cohort. (a) Placebo; and (b) levosimendan.

FIGURE 7.

Distribution of APACHE II score by treatment: the HARP-2 trial17 cohort. (a) Placebo; and (b) simvastatin.

As expected, mortality, in general, increased with increasing baseline risk measures for all trials. Figures 8–13 show the relationship between mortality and baseline risk measures in the VANISH trial,15 both for 28-day and for hospital mortality. Figures 14–19 show the same associations in the LeoPARDS trial16 and Figure 20 shows these associations for the HARP-2 trial. 17

FIGURE 8.

Mortality at day 28 (proportion) by APACHE II score: the VANISH trial15 cohort.

FIGURE 9.

Hospital mortality (proportion) by APACHE II score: the VANISH trial15 cohort.

FIGURE 10.

Mortality at 28 days (proportion) by APS-APII: the VANISH trial15 cohort.

FIGURE 11.

Hospital mortality (proportion) by APS-APII: the VANISH trial15 cohort.

FIGURE 12.

Mortality at day 28 (proportion) by R_calc. : the VANISH trial15 cohort.

FIGURE 13.

Hospital mortality (proportion) by R_calc. : the VANISH trial15 cohort.

FIGURE 14.

Mortality at day 28 (proportion) by APACHE II score: the LeoPARDS trial16 cohort.

FIGURE 15.

Hospital mortality (proportion) by APACHE II score: the LeoPARDS trial16 cohort.

FIGURE 16.

Mortality at 28 days (proportion) by APS-APII: the LeoPARDS trial16 cohort.

FIGURE 17.

Hospital mortality by APS-APII: the LeoPARDS trial16 cohort.

FIGURE 18.

Mortality at day 28 (proportion) by R_calc. : the LeoPARDS trial16 cohort.

FIGURE 19.

Hospital mortality (proportion) by R_calc. : the LeoPARDS trial16 cohort.

FIGURE 20.

Mortality at 28 days (proportion) by APACHE II score: the HARP-2 trial17 cohort.

Modified APACHE II risk of death model recalibrated

The predictive performance of APACHE II and R are shown in Table 3. After correction for overoptimism, model M3 without including the number of organ dysfunctions yielded the highest AUROC and DS (Table 4). The estimated parameters from the model are given in Table 5. In the original equation for R, a unit increase in the APACHE II score was associated with a 16% increase in the odds of hospital mortality, with scores from each variable having the same contribution to the prediction. In comparison, in model M3, the effect of unit increases in the APACHE II score on the odds of mortality ranged from a 17% decrease (temperature) to a 76% increase (pH). The baseline odds (i.e. the odds of mortality for a patient with an APACHE II score of zero and diagnostic category weight who was not admitted following emergency surgery) was lower for model M3 than in the original model for R (0.008 vs. 0.03). Therefore, the same OR would produce a much smaller absolute difference in model M3 than in the original model for R.

The VANISH trial heterogeneity of treatment effect assessment

In the primary analysis with APACHE II score as baseline risk of death measure, there was no evidence of HTE for vasopressin in either absolute terms [low-APACHE II group, RD 0.02 (95% CI –0.09 to 0.13); high-APACHE II group, RD 0.05 (95% CI –0.08 to 0.19); difference in RD 0.04 (95% CI –0.14 to 0.21)] or relative terms [low-APACHE II group, RR 1.09 (95% CI 0.64 to 1.86); high-APACHE II, group RR 1.15 (95% CI 0.80 to 1.64); ratio of RR 1.05 (95% CI 0.55 to 2.00)] (Figure 21). In the case of the secondary risk measures, the estimates of HTE for vasopressin were larger with wider CI for APS-APII (Figure 22) and smaller in magnitude for R (Figure 23).

FIGURE 21.

Forest plots for the RD and RR comparing 28-day mortality in treatment and control, by trial and APACHE II score subgroup. Reproduced from Santhakumaran and colleagues. 44 This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated. The figure includes minor additions and formatting changes to the original figure.

Reproduced from Santhakumaran and colleagues44 This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated. The figure includes minor additions and formatting changes to the original figure.

FIGURE 22.

Forest plots for the RD and RR comparing 28-day mortality in treatment and control, by trial and APS-APII subgroup. Reproduced from Santhakumaran and colleagues. 44 This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated. The figure includes minor additions and formatting changes to the original figure.

Reproduced from Santhakumaran and colleagues44 This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated. The figure includes minor additions and formatting changes to the original figure.

FIGURE 23.

Forest plots for the RD and RR comparing 28-day mortality in treatment and control, by trial and R subgroup. Reproduced from Santhakumaran and colleagues. 44 This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated. The figure includes minor additions and formatting changes to the original figure.

Reproduced from Santhakumaran and colleagues44 This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated. The figure includes minor additions and formatting changes to the original figure.

In the primary analysis with APACHE II score as baseline risk of death measure, there was no evidence of HTE for hydrocortisone in either absolute terms [low-APACHE II group, RD 0.02 (95% CI –0.12 to 0.17); high-APACHE II group, RD 0.06 (95% CI –0.10 to 0.21); difference in RD 0.03 (95% CI –0.18 to 0.25)] or in relative terms [low-APACHE II group, RR 1.11 (95% CI 0.62 to 1.99); high-APACHE II group RR 1.15 (95% CI 0.79 to 1.67); ratio of RR 1.04 (95% CI 0.52 to 2.08)]. In the case of the secondary risk measures, the estimates of HTE for hydrocortisone was similar for APS-APII (see Figure 22) and larger in magnitude for R (see Figure 23). Figures 21–23 were previously published by the authors in the paper by Santhakumaran and colleagues. 44

Heterogeneity of treatment effect was not observed when R was recalibrated either with controls only (Figure 24) or with the whole cohort (Figure 25), although subgroup differences were in the opposite direction for hydrocortisone.

FIGURE 24.

Forest plots for the RD and RR comparing 28-day mortality in treatment and control, by trial and R_recal. (recalibration model with only controls) subgroup.

FIGURE 25.

Forest plots for the RD and RR comparing 28-day mortality in treatment and control, by trial and R_recal. (recalibration model with whole cohort) subgroup.

The LeoPARDS trial heterogeneity of treatment effect assessment

For the primary analysis with APACHE II score as baseline risk of death measure there was no evidence of HTE for levosimenden in either absolute terms [low-APACHE II group, RD 0.05 (95% CI –0.04 to 0.15); high-APACHE II group, RD 0.04 (95% CI –0.08 to 0.16); difference in RD –0.02 (95% CI –0.17 to 0.14)] or in relative terms [low-APACHE II group, RR 1.34 (95% CI 0.78 to 2.31); high-APACHE II group, RR 1.09 (95% CI 0.84 to 1.41); ratio of RR 0.81 (95% CI 0.44 to 1.48)] (see Figure 21). For the secondary risk measures, the estimates of HTE for levosimenden were larger for APS-APII (see Figure 22) and in the opposite direction for R (see Figure 23).

Heterogeneity of treatment effect was not observed when R was recalibrated either with controls only (see Figure 24) or with the whole cohort (see Figure 25), although subgroup differences were in the opposite direction.

The HARP-2 trial heterogeneity of treatment effect assessment

For the primary analysis with APACHE II score as baseline risk of death measure, we observed HTE for simvastatin in absolute terms [low-APACHE II group, RD –0.15 (95% CI –0.22 to –0.07); high-APACHE II group, RD 0.19 (95% CI –0.01 to 0.39); difference in RD 0.34 (95% CI 0.12 to 0.55) (p = 0.02)] and in relative terms [low-APACHE II group, RR 0.45 (95% CI 0.28 to 0.72); high-APACHE II group, RR 1.61 (95% CI 0.95 to 2.71), ratio of RR 3.57 (95% CI 1.77 to 7.17)]. Simvastatin reduced mortality in the low-APACHE II group and increased mortality in the high-APACHE II group (see Figure 21). As raw data APACHE II score data were not available, we have not reported any secondary risk measures for the HARP-2 trial. 17

Serious adverse events and baseline risk

We plotted the proportions of serious adverse events in the low- and high-APACHE II groups in each trial to explore whether or not the pattern of adverse event distribution could explain any HTE in mortality. In all three RCTs, both in the intervention and controls trial arms, there was no pattern in serious adverse events that could explain HTE in mortality (Figure 26).

FIGURE 26.

Forest plots for the RD and RR comparing related serious adverse events in treatment and control, by trial and APACHE II score subgroup.

Heterogeneity of treatment effect assessment on continuous scale using regression

Differences were also smaller when HTE was assessed across the continuous range of APACHE II score [ratio of OR for 5-point increase in APACHE II 1.33 (95% CI 0.93 to 1.90)] (Table 6 and Figure 27).

FIGURE 27.

Heterogeneity of treatment effect assessment for APACHE II score as a continuous variable. Figures shows the estimated treatment effect with 95% CI bands from regression models for 28-day mortality including a treatment × APACHE II score interaction with (a) the VANISH trial15 (vasopressin); (b) the VANISH trial15 (hydrocortisone); (c) the LeoPARDS trial;16 and (d) the HARP-2 trial. 17

Sensitivity analyses

The results from sensitivity analyses were consistent with those from the main analyses for the VANISH trial15 and the LeoPARDS trial16 (see Table 6 and Figure 28). HTE was attenuated in the sensitivity analyses for the HARP-2 trial17 under different assumptions for the missing data [e.g. ratio of RR was 2.86 (95% CI 1.47 to 5.57) when we assumed that patients with missing APACHE II data were more likely to be high risk; all other results were less attenuated] (Table 7). Differences were also smaller when hospital mortality was used as the outcome [difference in RD 0.25 (95% CI 0.03 to 0.48), ratio of RR 2.34 (95% CI 1.31 to 4.18)] (see Figure 28).

FIGURE 28.

Forest plots for the RD and RR comparing hospital mortality in treatment and control, by trial and APACHE II score subgroup.

Exploratory analysis

Biomarker data (at least one biomarker at baseline) were available for 176 of 409 patients in the VANISH trial15 and 493 of 516 patients in the LeoPARDS trial. 16 Clinical characteristics at baseline are shown in Table 8. A summary of the biomarker data for both trials is shown in Table 9 (the VANISH trial15) and Table 10 (the LeoPARDS trial16), including details of values outside the limits of detection. As the limits varied by assay run, the mean limits for each biomarker are given.

Latent class analysis: the VANISH trial

The latent class modelling was carried out in three stages. The first stage includes only indicator variables (the biomarkers) and assumes a common variance across classes and zero covariance between indicators within class. In the second stage, clinical and demographic characteristics are added as covariates. In the third stage, the variance assumptions are relaxed. We first present the model results and then compare the results and model fit across all the models. Finally, we compare detailed results for a selection of candidate models. For each model we present the estimated distribution of the latent classes and the estimated class means for each indicator. Important indicators are those that have good separation (high between-class variability). As a measure of separation, we present the variance of the estimated class means. All indicators have been log-transformed and standardised to have a mean of zero and a SD of 1.

Comparing models derived from stages 1–3

The log-likelihood, class distributions, entropy, mean class probability, AIC and BIC are given in Table 11. Figure 29 shows how the log-likelihood, AIC and BIC change with the number of classes for each model stage. Across all models a similar set of indicators contributed to defining the classes. The estimated class means and class size were similar for the two-class models. Based on the AIC and BIC, the two-class stage 3b model (including covariates and allowing indicators to be correlated for individuals in the same class) appears to offer the best fit. We examined this model further, along with the three-class stage 3b model, which has a similar fit but an additional class, and the more parsimonious stage 2 models (two and three classes).

TABLE 11 - Model fit statistics for all LCA models in the VANISH trial15

Class statistics not applicable for models with only one class.

Model did not converge.

Stage	Number of classes	Log-likelihood	Class distribution		Entropy	Mean class probabilitya	AIC	BIC
			Estimated	Observeda
1	1a	–3937					7942	8049
	2	–3732	49/51	49/51	0.87	0.96/0.97	7568	7733
	3	–3673	29/44/27	29/45/26	0.82	0.93/0.89/0.94	7486	7708
	4	–3638	21/13/34/32	22/12/35/31	0.84	0.91/0.90/0.88/0.95	7451	7730
2	1a	–3806					7714	7874
	2	–3596	51/50	49/51	0.88	0.97/0.96	7378	7671
	3	–3514	37/28/35	36/26/38	0.86	0.95/0.98/0.90	7298	7722
	4	–3465	36/17/24/23	35/17/23/24	0.88	0.95/0.95/0.93/0.91	7281	7838
3a	1a	–3806					7714	7875
	2	–3573	51/49	49/51	0.89	0.97/0.96	7365	7711
	3	–3476	49/29/23	47/28/24	0.93	0.98/0.96/0.95	7290	7822
	4b
3b	1a	–3621					7453	7787
	2	–3454	46/54	44/56	0.85	0.96/0.95	7025	7671
	3	–3375	32/46/22	31/45/24	0.89	0.95/0.94/0.98	7130	7729
	4	–3330	25/38/23/15	24/37/23/16	0.91	0.96/0.95/0.95/0.96	7125	7855
3c	1a	–3621					7453	7787
	2	–3393	56/44	55/45	0.90	0.97/0.96	7226	7918
	3	–3307	48/23/29	48/22/31	0.93	0.97/0.99/0.97	7282	8333
	4b

FIGURE 29.

Plots of model fit indicators in the VANISH trial:15 (a) log-likelihood; (b) AIC; and (c) BIC.

Figure 30 shows the distribution of each indicator by class for each of the candidate models, assigning individuals to their modal class (i.e. the class for which they had the highest posterior class probability). The indicators are ordered by the p-value from a Kruskal–Wallis non-parametric test comparing the distribution across the classes. The test was performed separately for each model and the average p-value taken to get an approximate ordering of importance of the indicators. The results are consistent across the models, with differences observed only for indicators with less separation across the classes. For example, NT-proBNP is highest in class 2 for the stage 2 three-class model, but highest in class 3 for the stage 3 three-class model. Figures 41–43 show separation plots for each of the candidate models reported in stages 1–3. These show the distribution of each of the indicators, with one line representing an individual, coloured according to the modal class. They show how well separated the classes are for each indicator. Only the top few indicators are clearly separated. For the three-class models, stage 3b models appear slightly better separated than their stage 2 counterpart for the most important indicators. There is little difference between the two-class models.

FIGURE 30.

Box plot showing distribution of indicators by class for candidate models, assigning to modal class, in the VANISH trial. 15 The box shows the 25th, 50th and 75th percentiles, the whiskers show the extent of data within 1.5 × IQR and the dots show data points beyond these limits. P/F, PaO₂/FiO₂.

In summary, these comparisons suggest that the simpler two-class model from stage 2 is appropriate, given the minimal differences in indicator distribution, separation and class assignment. The more complex three-class models create a ‘middle’ class that is a mixture of the other two classes, and so does not give a substantively different interpretation. Therefore, the two-class model from stage 2 is used for all our subsequent analysis.

Analysis by modal class: the VANISH trial

There was an almost even split between classes (90 individuals assigned to class 1 and 86 individuals assigned to class 2). The clinical characteristics by latent class in the VANISH trial15 are shown in Table 12, biomarker values are shown in Table 13 and clinical outcomes by class are shown in Table 14. The classes in the final model will be referred to as subphenotype 1 (i.e. class 1) and subphenotype 2 (i.e. class 2) from hereon in the manuscript when referring to the VANISH trial. 15

TABLE 12 - Baseline characteristics by assigned class in the LeoPARDS trial16 and the VANISH trial15

Characteristic	LeoPARDS trial16 (N = 493)				VANISH trial15 (N = 176)
	Class 1 (n = 191)	Class 2 (n = 247)	Class 3 (n = 55)	Missing (n)	Class 1 (n = 90)	Class 2 (n = 86)	Missing (n)
Age (years), median (IQR)	68 (57–77)	69 (62–76)	65 (51–73)	0	65.5 (54–77)	64.5 (53–76)	0
Male, n (%)	108 (56.5)	138 (55.9)	28 (50.9)	0	58 (64.4)	54 (62.8)	0
Ethnicity, n (%)				0			0
White	180 (94.2)	229 (92.7)	52 (94.6)		77 (85.6)	69 (80.2)
Black	3 (1.6)	5 (2)	2 (3.6)		8 (8.9)	6 (7)
Asian	7 (3.7)	11 (4.5)	1 (1.8)		4 (4.4)	9 (10.5)
Other	1 (0.5)	2 (0.8)	0 (0)		1 (1.1)	2 (2.3)
BMI (kg/m2), median (IQR)	26.7 (23–30.8)	27.3 (23.4–30.7)	27.8 (24.2–33.6)	6	24.7 (22.2–31.6)	26.2 (22.6–31.1)	9
Comorbidities, n (%)
NYHA IV	1 (0.5)	4 (1.6)	0 (0)	0	0 (0)	0 (0)	0
Severe COPD	10 (5.2)	11 (4.5)	2 (3.6)	0	7 (7.8)	3 (3.5)	0
Chronic renal failure	18 (9.4)	15 (6.1)	2 (3.6)	0	2 (2.2)	6 (7)	0
Cirrhosis	6 (3.1)	3 (1.2)	0 (0)	0	4 (4.4)	7 (8.1)	0
Immunocompromised	14 (7.3)	22 (8.9)	9 (16.4)	0	2 (2.2)	9 (10.5)	0
Site of infection, n (%)				3			1
Lung	106 (55.5)	74 (30.1)	12 (21.8)		42 (48.3)	32 (37.2)
Abdomen	45 (23.6)	111 (45.1)	25 (45.5)		17 (19.5)	18 (20.9)
Urine	13 (6.8)	13 (5.3)	3 (5.5)		10 (11.5)	18 (20.9)
Primary bacteraemia	0 (0)	6 (2.4)	4 (7.3)		2 (2.3)	1 (1.2)
Neurological	3 (1.6)	2 (0.8)	0 (0)		4 (4.6)	0 (0)
Soft tissue or line	10 (5.2)	10 (4.1)	6 (10.9)		3 (3.5)	(3.5)
Other	14 (7.3)	30 (12.2)	5 (9.1)		9 (10.3)	14 (16.3)
SOFA score (points), median (IQR)	7 (6–8)	8 (7–10)	9 (7–11)	22	6 (4–8)	8 (5–10)	16
APACHE II score (points), median (IQR)	24 (21–30)	26 (21–31)	27 (22–30)	1	24 (18–29)	23.5 (20–30)	0
Post-surgical admission, n (%)	49 (25.7)	114 (46.2)	17 (30.9)	0	17 (18.9)	9 (10.5)	0

TABLE 13 - Biomarker data by class, with study participants assigned by highest posterior class probability

Parameter	LeoPARDS trial16 (N = 493)				VANISH trial15 (N = 176)
	Class 1, median (IQR)	Class 2, median (IQR)	Class 3, median (IQR)	n (%) outside limits	Class 1, median (IQR)	Class 2, median (IQR)	n (%) outside limits
Organ dysfunction
PaO2/FiO2 ratio (kPa)	29.1 (22–39.7)	29.3 (20.2–39.4)	25.8 (16.6–36)	0/491 (0)	32.5 (21–43.5)	21.1 (14.6–34.9)	0/169 (0)
Creatinine (µmol/l)	107 (69–166)	151 (107–231)	173 (137–295)	0/491 (0)	91.5 (67–163)	140 (106–270)	0/176 (0)
Platelets (× 109/l)	243 (182–350)	203 (131–294)	136 (76–215)	0/490 (0)	206 (145–335)	150 (83–246)	0/171 (0)
Bilirubin (µmol/l)	12 (7–19)	17 (10–30)	17 (9–31)	0/483 (0)	12 (7–23)	16.5 (11–42)	0/156 (0)
Inflammation markers (pg/ml)
IL-1β	0.915 (0.651–1.41)	1.53 (0.948–2.88)	7.96 (2.93–11.3)	43/486 (8.8)	4.9 (4.34–6.25)	11.2 (7.71–19.4)	109/162 (67.3)
IL-6	232 (92–481)	1588 (583–3874)	19,582 (11,926–27,584)	34/490 (6.9)	426 (175–1376)	6385 (2277–19,641)	15/162 (9.3)
IL-8	48.2 (30.4–84.8)	257 (159–516)	3015 (1252–7336)	4/490 (0.8)	64.8 (28.3–173)	1075 (225–3293)	16/162 (9.9)
IL-10	26.3 (17.4–49.8)	123 (66.1–205)	554 (314–1429)	1/490 (0.2)	15.1 (8.27–32.9)	159 (66–446)	16/162 (9.9)
IL-17	6.61 (4.71–10.1)	9.79 (6.52–19.5)	21.5 (8.1–49.9)	9/486 (1.9)	15.6 (7.36–22.6)	18 (10.6–38.8)	108/162 (66.7)
IL-18	559 (373–996)	804 (565–1278)	1065 (724–1759)	16/486 (3.3)	434 (163–595)	562 (331–836)	3/162
Leucocytes (pg/ml)
Myeloperoxidase	332,192 (204,356–581,390)	489,569 (291,316–987,773)	541,091 (340,269–1,405,943)	47/486 (9.7)	264,802 (122,264–433,825)	675,769 (448,047–1,195,365)	3/168
sICAM	271,181 (168,034–414,706)	311,864 (187,557–515,209)	432,439 (298,658–886,124)	23/486 (4.7)	243,841 (165,654–326,461)	341,912 (208,505–550,529)	6/168
Endothelial injury (pg/ml)
ANG II	3197 (1906–5419)	7487 (4491–14,867)	13,040 (7373–23,168)	15/486 (3.1)	3162 (1503–5551)	6592 (3681–11,564)	3/168 (1.8)
Cardiovascular
Troponin (ng/l)	62 (16.8–536)	77.8 (23.7–381)	139 (45.2–589)	0/483	21 (6–94)	149 (31–729)	0/95
NT-proBNP (pg/ml)	9054 (3318–17,410)	10,269 (4922–23,317)	18,406 (9844–31,718)	36/492 (7.3)	3611 (1238–7416)	8130 (3590–16,890)	3/168 (1.8)
Other markers
sTNFR1 (pg/ml)	5939 (3923–9802)	13,457 (8749–20,337)	18,099 (12,379–27,759)	0/492	3856 (2064–5555)	8315 (5743–12,645)	2/168 (1.2)
Lactate (mmol/l)	1.5 (1–2.1)	2.6 (1.8–4)	5.2 (3–7)	0/490	1.8 (1.2–2.5)	3.5 (2.3–5.3)	0/172
CCL2 (pg/ml)	384 (272–592)	995 (676–1590)	4049 (3105–5621)	6/490 (1.2)

TABLE 14 - Clinical outcomes by class in the LeoPARDS trial16 and the VANISH trial15

Chi-squared test.

Kruskal–Wallis test.

In patients not in renal failure at baseline.

In patients who die or experience some renal failure by day 28.

Mann–Whitney U-test.

Trial outcome	Class			p-value for difference
	1	2	3
LeoPARDS16
3-month survival, n/N (%)	132/189 (69.8)	155/246 (63.0)	23/55 (41.8)	0.001a
28-day survival, n/N (%)	143/190 (75.3)	165/247 (66.8)	26/55 (47.3)	< 0.001a
Mean daily SOFA score (points), mean (SD)	4.93 (2.88)	6.67 (3.94)	9.97 (4.60)	< 0.001b
VANISH15
28-day renal failure-free survival, n/N (%)c	51/75 (68.0)	32/60 (53.3)		0.08a
28-day survival, n/N (%)	72/90 (80.0)	57/86 (66.3)		0.04a
Renal failure-free days, median (IQR)d	19 (3–26)	8 (0–23)		0.03e

There was no evidence that treatment effects varied by subphenotype for any of the outcomes of the VANISH trial15 (Figures 31 and 32). The effect of vasopressin on renal failure-free survival at 28 days compared with noradrenaline was in opposite directions in class 1 [i.e. 10% reduction in survival (95% CI –31% to 11%)] compared with subphenotype 2 [i.e. 10% increase in survival (95% CI –16% to 35%)], but the CI for the subgroup difference was wide [difference in RD 20% (95% CI –13% to 53%)]. Point estimates for the treatment effect showed a consistent direction of subgroup differences (i.e. all RDs were positive, indicating that treatments were more likely to benefit participants in class 2). This was also seen for renal failure-free days. For subphenotype 1, the median in the vasopressin group was 10 days lower (95% CI –23 to 3 days) than in the noradrenaline group, whereas for class 2 those in the vasopressin group had a median of 6 more renal failure-free days (95% CI –8 to 20 days). The test for subphenotype–treatment interactions was not statistically significant. The sensitivity analysis gave very similar results, with 97% agreement in subphenotype assignment for the VANISH trial15 (see Table 24).

FIGURE 31.

No evidence for treatment effect variation by class for any of the survival outcomes of the VANISH trial15 and the LeoPARDS trial. 16 RF, renal failure.

FIGURE 32.

No evidence for treatment effect variation by class for any of the organ dysfunction outcomes of the VANISH trial15 and the LeoPARDS trial. 16

Latent class analysis: the LeoPARDS trial

Comparing models derived from stages 1–3

A similar set of indicators were important for defining the latent classes, with some overlap with those found in the VANISH trial. 15 The log-likelihood, class distributions, entropy, mean class probability, AIC and BIC are given in Table 15. We present plots showing how the log-likelihood, AIC and BIC change with the number of classes for each model stage in Figure 33.

TABLE 15 - Model fit statistics for all models: the LeoPARDS trial16

Class statistics not applicable for models with only one class.

The five-class model did not converge in stages 2–3c.

Stage	Number of classes	Log-likelihood	Estimated class distribution	Observed distributiona	Entropy	Mean class probabilitya	AIC	BIC
1	1a	–13,058					26,188	26,339
	2	–12,285	61/39	61/39	0.89	0.98/0.96	24,681	24,912
	3	–11,942	33/14/53	32/14/54	0.89	0.95/0.97/0.95	24,031	24,342
	4	–11,800	35/23/13/29	35/23/13/30	0.86	0.95/0.89/0.98/0.88	23,787	24,178
	5b	–11,676	22/30/17/10/20	22/30/17/10/20	0.86	0.93/0.89/0.91/0.98/0.90	23,576	24,047
2	1a	–12,657					25,423	25,648
	2	–11,910	58/42	57/43	0.89	0.98/0.96	24,020	24,438
	3	–11,555	32/12/56	31/12/58	0.91	0.95/0.99/0.95	23,402	24,013
	4	–11,373	26/12/21/41	25/11/20/43	0.89	0.96/0.98/0.93/0.92	23,129	23,932
3a	1a	–12,657					25,423	25,648
	2	–11,770	54/46	53/47	0.88	0.96/0.97	23,776	24,270
	3	–11,413	30/51/19	29/51/20	0.90	0.97/0.95/0.97	23,189	23,950
	4	–11,174	31/29/16/24	31/29/15/25	0.90	0.96/0.93/0.98/0.94	22,840	23,868
3b	1a	–12,470					25,091	25,404
	2	–11,796	58/42	57/43	0.87	0.97/0.96	23,833	24,339
	3	–11,452	39/11/50	39/11/50	0.91	0.94/0.98/0.96	23,238	23,936
	4	–11,285	31/11/33/25	31/11/32/27	0.88	0.95/0.98/0.93/0.90	22,997	23,887
3c	1a	–12,470					25,091	25,404
	2	–11,641	54/46	53/47	0.88	0.96/0.97	23,602	24,271
	3	–11,283	37/43/20	36/43/21	0.90	0.97/0.94/0.98	23,055	24,079
	4	–11,004	33/30/16/21	32/30/16/23	0.91	0.97/0.93/0.98/0.95	22,749	24,128

FIGURE 33.

Plots of model fit indicators in the LeoPARDS trial:16 (a) log-likelihood; (b) AIC; and (c) BIC.

With the exception of the stage 3c model, the information criteria continued to decrease as the number of classes increased. There was negligible improvement after three classes and the BIC for the stage 3c model increased. Fit statistics were close for stages 2–3c and so we examined the three-class model for all of these stages in more detail.

Figure 34 shows the distribution of each indicator by class for each of the candidate models, assigning individuals to their modal class (i.e. the class for which they had the highest posterior class probability). The results are very similar across the models. The most notable difference is for values of IL-6 in class 3, which are more dispersed in the models and allow for the variance to differ across the classes (i.e. stages 3a and 3c).

FIGURE 34.

Box plot showing distribution of indicators by class for candidate models, assigning to modal class, in the LeoPARDS trial. 16 The box shows the 25th, 50th and 75th percentiles, the whiskers show the extent of data within 1.5 × IQR and the dots show data points beyond these limits. P/F, PaO₂/FiO₂.

Figures 44–47 show separation plots for each of the candidate models. These show the distribution of each indicators, with one line representing an individual, coloured in accordance with the modal class. They show how well separated the classes are for each indicator. For stages 3a and 3c there are a few low biomarker values that belong to subjects in class 3.

In summary, there are minimal difference in indicator distribution. Models 3a and 3c show some indication of poor separation, and model 2 may be of insufficient complexity, given what is known about the correlations between the indicators. Therefore, model 3b was selected.

Analysis by modal class: the LeoPARDS trial

The baseline clinical and demographic characteristics by classes are shown in Table 12. The final model assigned 191 individuals to class 1, 247 individuals to class 2 and 55 individuals to class 3. As standardised values can be difficult to interpret, the median and IQR for each indicator on the original scale are shown in Table 13 for individuals assigned to each class, based on posterior class probability. The classes in the final model will be referred to as subphenotype 1 (i.e. class 1), subphenotype 2 (i.e. class 2) and subphenotype 3 (i.e. class 3) from hereon in the manuscript.

Survival varied by subphenotype (p = 0.001). In particular, survival was lower in subphenotype 3 (23/55, 41.8%) than in the other subphenotypes [subphenotype 1, 132/189 (69.8%); subphenotype 2, 155/246 (63.0%)]. Similar results were seen for survival to 28 days (see Table 14). The mean daily SOFA score also increased with subphenotype, with score almost twice as high in subphenotype 3 than it was in subphenotype 1.

There was no evidence that treatment effects varied by subphenotype for any of the outcomes of the LeoPARDS trial. 16 Survival was lower in the levosimendan group for all classes (although the difference was not statistically significant in any subphenotype), with no apparent trend across the subphenotypes (see Figure 31). Mean daily SOFA score was higher in the levosimendan group than in the placebo group in subphenotype 1 (RD 0.84, 95% CI 0.02 to 1.65) and in subphenotype 2 (RD 0.88, 95% CI –0.08 to 1.83) (see Figure 32). There was no evidence of treatment difference in subphenotype 3 (RD –0.14, 95% CI –2.55 to 2.28). However, the differences in treatment effect comparing classes were not statistically significant.

A multinomial logit model with IL-6, IL-8, IL-10 and CCL2 as predictors gave a sensitivity of around 0.9 and a specificity of ≥ 0.9 for all subphenotypes (Figure 35). In particular, for subphenotype 3, which is perhaps of most interest for the purposes of identifying a trial population, the specificity was 0.98. The model coefficients are shown in Table 16. The sensitivity analysis gave very similar results, with 94% agreement in subphenotype assignment for the LeoPARDS trial16 (see Table 24).

FIGURE 35.

Class-specific (a) sensitivity, (b) specificity and (c) c-statistics for multinomial logit models with increasing number of predictors in the LeoPARDS trial. 16

TABLE 16 - Model coefficients from the multinomial regression model in the LeoPARDS trial16

Predictor	Class 2 vs. class 1	SE	Class 3 vs. class 1	SE
	Log-OR for a 1 SD increase		Log-OR for a 1 SD increase
IL-6	1.404	0.356	3.803	0.771
IL-8	2.477	0.472	4.09	0.735
IL-10	1.942	0.346	2.514	0.639
CCL2	1.314	0.328	3.976	0.853
Baseline odds	2.236	0.301	–8.264	1.952

Latent class analysis: the HARP-2 trial

Parts of this section, which presents data on ARDS subphenotypes from the HARP-2 trial,17 includes information based on our previous publication by Calfee and colleagues. 31

Comparison of phenotypic features and outcomes between subphenotypes

Subphenotype 2 had clinical and biological features similar to those found in our prior studies and consistent with a hyperinflammatory phenotype. Specifically, when compared with subphenotype 1, patients in subphenotype 2 had higher values of sTNFR1 and IL-6, lower platelet counts (Figure 36) and more vasopressor use (p < 0.001). Age and sex were similar across the subphenotypes. Although the distribution of direct and indirect ARDS risk factors was significantly different across the two subphenotypes (p < 0.0001), the most common ARDS risk factors of sepsis, pneumonia and aspiration were highly prevalent among both groups, as in our prior work. 23 In addition, subphenotype 2 patients had fewer ventilator-free days (median of 2 vs. 18 days; p < 0.0001), fewer non-pulmonary organ failure-free days (median of 15 vs. 27 days; p < 0.0001) and higher 28-day mortality (39% vs. 17%; p < 0.0001) than subphenotype 1 patients (Table 18).

FIGURE 36.

Differences in standardised values of each continuous variable by subphenotype in the HARP-2 trial. 17 Standardised means by class. Variables with maximum positive separation are shown on the left (i.e. hyperinflammatory subphenotype higher than hypoinflammatory subphenotype) and variables with maximum negative separation on the right (i.e. hyperinflammatory subphenotype lower than hypoinflammatory subphenotype). The y-axis represents standardised variable values. A value of +1 for the standardised variable signifies that the mean value for a given subphenotype was 1 SD higher than the mean value in the cohort as a whole. The mean values are joined by lines to facilitate display of subphenotype profiles. P/F, PaO₂/FiO₂.

TABLE 18 - Clinical outcomes by subphenotype in the HARP-2 trial17

Outcome	Class		p-value
	1 (n = 354)	2 (n = 186)
28-day mortality, n (%)	59 (17)	73 (3)	< 0.0001
90-day mortality, n (%)	78 (22)	87 (46)	< 0.0001
Ventilator-free days, median (25–75%)	2 (0–17)	18 (0–23)	< 0.0001
Non-pulmonary organ failure-free days, median (25–75%)	27 (21–28)	15 (0–25)	< 0.0001

Survival benefit observed with simvastatin in subphenotype 2

The original trial found no difference in 28-day survival curves between placebo and simvastatin (p = 0.20). The subphenotype 2 had a better 28-day survival (p < 0.0001) (Figure 37) and better 90-day survival (p < 0.0001) (Figure 38) for overall comparison (p = 0.03 for subphenotype 2 simvastatin vs. placebo and p = 0.21 with Bonferroni correction).

FIGURE 37.

Kaplan–Meier survival curves. (a) 28-day patient survival; and (b) 90-day patient survival in HARP-2 trial,17 stratified by subphenotype and treatment (simvastatin vs. placebo).

FIGURE 38.

Survival curves stratified by mean APACHE II score. Product-limit survival estimates. (a) 28-day survival; and (b) 90-day patient survival in the HARP-2 trial,17 stratified by ARDS subphenotype and treatment (simvastatin vs. placebo).

In contrast to the curves stratified by subphenotype and treatment, survival curves stratified by ARDS severity (i.e. PaO₂/FiO₂ ratio) and treatment were not significantly different (p = 0.12). Survival curves stratified by APACHE II score (dichotomised at the median) and treatment revealed differences in survival by APACHE II score, but no differential effect of treatment in either the high- or low-APACHE II group (see Figure 38).

Mortality at 28 days was 13% lower in subphenotype 2 patients treated with simvastatin than in subphenotype 2 patients treated with placebo (32% vs. 45%). In contrast, 28-day mortality was similar in patients in subphenotype 1 regardless of treatment assignment (16% vs. 17%). The interaction between treatment and subphenotype for mortality was not statistically significant (p = 0.14).

In the original trial, time to unassisted breathing did not differ significantly between simvastatin- and placebo-treated patients, although a trend favouring simvastatin was observed (hazard ratio 0.84; p = 0.09). When stratified by subphenotype and treatment, time to unassisted breathing differed significantly (p < 0.0001) (Figure 39). However, the difference in the curves between subphenotype 2 patients treated with simvastatin and placebo was not statistically significant (p = 0.10). Among subphenotype 2 patients, median ventilator-free days were numerically higher in the simvastatin-treated patients than in placebo-treated patients (7 days vs. 0 days). This is in contrast to patients in subphenotype 1, among whom the median number of ventilator-free days was the same regardless of treatment (18 days in each). However, the interaction between treatment and subphenotype in regression models was not statistically significant (p = 0.15).

FIGURE 39.

Time to unassisted breathing over 28 days in the HARP-2 trial. 17 Data are stratified by subphenotype and treatment condition, from Fine–Gray competing risks model. 59 HR, hazard ratio.

Main findings

We assessed whether or not HTE could contribute to the results in three recent ICU RCTs, using four different multivariable baseline risk of death models, which included well-established risk factors for acute mortality for sepsis and ARDS as covariates. There was considerable within-trial variation in the baseline risk of death in all three RCTs. We did not find consistent evidence for HTE as an explanation for the original trial results in all three ICU RCTs we assessed, as tests for HTE were inconclusive. We observed that detection of HTE in RCTs may be influenced by the baseline risk model specification, as illustrated by differences in HTE effects seen in the LeoPARDS trial. 16 The lack of evidence of relative HTE may be in part due to low power to detect interaction effects, particularly for the HARP-2 trial,17 in which mortality was lower among those allocated to simvastatin in all risk quartiles except the highest. Parts of this section, which presents data on ARDS subphenotypes from the HARP-2 trial,17 includes information based on our previous publication by Calfee and colleagues. 31

Explanation of key findings

There are a number of possible reasons why we did not observe HTE consistently in our analyses. All three trials we assessed have many features of explanatory trials,60 which by their design limit HTE in comparison with pragmatic trials. Therefore, demonstrable HTE is less likely in these trials, although its evaluation remains important. Our findings may be true in that HTE may be less marked in sepsis and ARDS, in which many ‘minimal causes’ of mortality may be contributory,61 than in illnesses such as retroviral disease. 62 It could be that the effects of the treatments we assessed are small and of limited variability, resulting in minimal HTE. In other words, when HTE is assessed for an intervention using data from a single trial we are unlikely to detect it unless HTE effects are large. This generates an argument to assess HTE using trials of the similar treatment–condition combination or of the same condition, and a broader group of treatments with similar enough mechanism of treatment effect or to consider intervention-specific multivariable models.

Comparison with published literature

A key comparison to consider is the contrasting results with RCT simulations by Iwashyna and colleagues. 24 Their simulations assumed that the trial participants’ odds of 30-day mortality will be influenced by severity of acute respiratory failure, comorbid conditions, the extent to which the treatment reduces mortality from the primary illness and the treatment’s fatal adverse effect rates. Although we used 28-day mortality for our primary analysis, we considered baseline risk as a function of acute illness severity using the total APACHE II score. The data in our trials do not follow the patterns described by Iwashyna and colleagues,24 as we did not have constant relative treatment effects or constant harms, and mortality patterns differed from those predicted by their simulation model. Recently, Semler and colleagues63 reported a pragmatic, cluster-randomised, multiple-crossover trial of saline compared with balance crystalloids in critically ill patients, and found no difference in the primary outcome of major adverse kidney events within 30 days, but a positive result in the composite secondary outcomes of death from any cause, new renal replacement therapy and persistent renal dysfunction. The authors then reported presence of HTE for the primary trial outcome, when assessed with a multivariable model specifically calibrated for the primary outcome. 64 Therefore, it is plausible that the closer a trial is to showing a difference, the greater the chance of finding HTE with multivariable models, which might also explain our findings of a lack of HTE.

Strengths and weakness

We explored HTE in sepsis and ARDS for four different treatments and using four different multivariable models. The primary baseline risk measure, APACHE II score, is an established, validated predictor of mortality in this population. Three variations on this measure were investigated to check the consistency of the results, along with several sensitivity analyses. We used a composite risk score (APACHE II) for its superior performance for baseline risk estimation, as highlighted by Kent and colleagues18 and as recommended for future studies of HTE assessment. 24 Furthermore, lead time bias may influence the impact that APACHE II score might have on the outcome. This may be less relevant here, as all the data sets were from RCTs. None of the RCTs included in this study had 28-day mortality as the primary outcome. Therefore, it is possible that we were underpowered to detect HTE, if it existed. The primary outcomes were not suitable for HTE analysis because they were continuous rather than binary and without an appropriate baseline measure, although the existing HTE framework could be adapted for some continuous outcomes, such as change from baseline organ dysfunction. Similarly, there may have been insufficient numbers to demonstrate non-linear risk and HTE effects that we investigated as sensitivity analyses.

Main findings in the VANISH trial and the LeoPARDS trial

The VANISH trial15 two-class (subphenotype) model best represented the data and generated two sepsis subphenotypes with a split of 90 (subphenotype 1) and 86 (subphenotype 2) individuals. The subphenotype 2 individuals had greater inflammation, with higher concentrations of IL-1β, IL-6, IL-8, IL-10, MPO, ANG II, troponin, NT-proBNP and sTNFR1. The class 2 individuals also had reduced survival and fewer renal failure-free days. There were no treatment effect differences between the two classes for corticosteroids or vasopressin.

The LeoPARDS trial16 three-class (subphenotype) model best represented the data and generated three sepsis subphenotypes with a split of 191 individuals to class 1, 247 individuals to class 2 and 55 individuals to class 3. The subphenotype 3 individuals had greatest inflammation, with higher concentrations of IL-1β, IL-6, IL-8, IL-10, IL-17, ANG II, troponin, NT-proBNP, CCL2 and sTNFR1. The subphenotype 3 individuals also had the highest SOFA score and reduced survival at 90 days. There were no treatment effect differences between classes for levosimendan, albeit survival was lower in the levosimendan group for all classes. A multinomial logit model with IL-6, IL-8, IL-10 and CCL2 as predictors gave a sensitivity of around 0.9 and a specificity of ≥ 0.9 for all classes. The differences in number of classes between the two trials may be related to differences in eligibility criteria, differences in sample sizes and the timing of measurements.

Comparison with published literature

Sepsis phenotypes can be categorised into clinical phenotypes that are identified using clinical and commonly acquired laboratory data (such as leucocyte count and C-reactive protein) and into molecular phenotypes that are identified using leucocyte gene expression data. This approach to classify phenotypes into clinical phenotype and molecular phenotypes has been reported in asthma. 65,66 Studies to date have identified four37,38 to six39 sepsis clinical phenotypes and two21,67 to four20 sepsis molecular phenotypes.

Using clinical data from the PROWESS Shock RCT,68 Gårdlund and colleagues39 reported six different sepsis clinical phenotypes using LCA, which were (1) uncomplicated septic shock, (2) pneumonia with ARDS, (3) postoperative abdominal shock, (4) severe septic shock, (5) pneumonia with ARDS and multiple organ dysfunction syndrome and (6) late septic shock. They did not report any biomarker data in their LCA to make direct comparisons with our study. Knox and colleagues38 reported four phenotypes using self-organising maps and grouped them into (1) shock with elevated creatinine; (2) minimal multiorgan dysfunction syndrome; (3) shock with hypoxaemia and altered mental status; and (4) hepatic disease. Seymour and colleagues37 reported four similar sepsis clinical phenotypes using observational cohort data and referred to them as (1) alpha (i.e. fewest abnormal laboratory test results, least organ dysfunction and lowest mortality), (2) beta (i.e. older patients, high prevalence of chronic illnesses and kidney dysfunction), (3) gamma (greater inflammation and pulmonary dysfunction) and (4) delta (i.e. with liver dysfunction and septic shock). In addition, Seymour and colleagues37 correlated the phenotype data with biomarkers, validated these findings within RCT data and simulated for treatment effect heterogeneity. They neither did a direct interaction test with treatment effect nor used the biomarkers we used in their LCA models. Based on pan-leucocyte gene expression data, Davenport and colleagues21 present two sepsis molecular phenotypes named sepsis response signatures [i.e. sepsis response signature 1 (SRS1) and sepsis response signature 2 (SRS2)], with the SRS1 phenotype being more immunosuppressed and having greater mortality than SRS2. A two-subset model based on neutrophils gene expression was reported by Maslove and colleagues,67 with subset 1 showing evidence of greater inflammation and toll-like receptor signalling. Based on pan-leucocyte gene expression data, Scicluna and colleagues20 reported four sepsis molecular phenotypes and subphenotypes assessment, albeit they do not strictly meet the endotype definition. We cannot provide a direct comparison with these results as we do not use transcriptome data in our LCA. Recently, Antcliffe and colleagues69 tested whether or not SRS1 and SRS2 sepsis molecular phenotypes respond differently to corticosteroids, using the VANISH trial15 data. The authors show that the relatively more immunocompetent phenotype (i.e. SRS2) may be harmed from corticosteroid therapy. In contrast, we did not observe any treatment effect differences between the two sepsis subphenotypes identified when using cytokines and clinical data in the VANISH trial. 15

Strengths and weakness

This study has several strengths. First, the LCA plan was defined prior to knowing the results of biomarker measurements. Second, we used data from two high-quality sepsis trials. 15,16 Last, the subphenotypes assessment followed a systematic plan and the model that best explained the heterogeneity within the trial cohorts was used.

The limitations are that this is a post hoc analysis, albeit with an explicit peer-reviewed hypothesis and an analysis plan that was set up prior to measurement of any biomarker data. We report two subphenotypes from the VANISH trial15 and three from the LeoPARDS trial,16 which is partly driven by the sample size. Not all patients in both trials had biomarker data. Some of the biomarker measurements were outside the measurement range. In addition, because of the nature of latent class models, there is uncertainty as to the similarities between two subphenotypes identified in the VANISH trial15 and the three subphenotypes identified in the LeoPARDS trial,16 although the key biomarkers appear similar between subphenotype 2 in the VANISH trial15 and subphenotype 3 in the LeoPARDS trial. 16

Main findings in the HARP-2 trial

The HARP-2 trial17 analysis shows two novel findings. First, two distinct ARDS subphenotypes with features similar to those previously reported were identified for the first time in a non-US patient population. Importantly, this was achieved using a different and much smaller set of clinical and biomarker data than in previous studies. Second, and more importantly, these two subphenotypes of ARDS responded differently to randomly assigned simvastatin, with evidence of improved survival at both 28 and 90 days uniquely among patients with a ‘hyperinflammatory’ subphenotype of ARDS. The finding that patients with a hyperinflammatory ARDS subphenotype preferentially responded to randomly assigned simvastatin treatment has biological plausibility based on the presumed mechanism of action of statins in ARDS. Statins reduce lung inflammation and injury in both animal models of ARDS and pre-clinical human experimental studies,70 and also have endothelium-stabilising properties. Therefore, patients with a higher degree of systemic inflammation, such as those in the hyperinflammatory subphenotype, would seem to be most likely to respond to this therapy.

Comparison with published acute respiratory distress syndrome literature

To date, there have been five RCT cohort reanalyses, including HARP-2 trial analysis23,31,40–42 and one observational cohort study. 22 All RCT reanalyses were conducted using LCA. 36 All RCTs show two phenotypes: (1) a non-inflammatory subphenotype that accounts for two-thirds of the trial population and (2) a less prevalent hyperinflammatory subphenotype that accounts for the remaining one-third of the population. The hyperinflammatory phenotype is associated with higher concentrations of inflammatory biomarkers (i.e. IL-6, IL-8 and sTNFR1), more acidosis and greater prevalence of shock, as suggested by vasopressor requirements. 36 Mortality is significantly higher among patients with the hyperinflammatory subphenotype than among those with the non-inflammatory phenotype. The observational cohort study also highlights a two-subphenotype model (i.e. equal prevalence of reactive and uninflamed populations). 36 When leucocyte gene expression was compared between the reactive and uninflamed subphenotypes in a cohort of 210 patients with sepsis and ARDS, 128 patients had a reactive phenotype and 82 patients had an uninflamed phenotype, with significant differences in 3332 of the 11,443 (29%) transcripts between the phenotypes. These findings highlight the need for a prospective study to validate the two ARDS subphenotypes using key discriminant markers. 71 These findings indicate that these subphenotypes are consistent across geographical sites and are robust to variations in specific data collected, enhancing the generalisability of previous studies.

Strengths and weakness

The HARP-2 trial17 analysis is a post hoc analysis with attendant limitations, including lack of statistical power, as implied by the statistical test for interaction in the analyses of 28-day mortality (p = 0.14). 72 We acknowledge that this analysis was not planned as part of the original trial design because the trial was designed before the descriptions of ARDS subphenotypes. 73 In addition, because of the nature of latent class models, it is not possible to prove that the two subphenotypes identified in the HARP-2 trial17 are ‘the same’ as the two subphenotypes identified in previous studies. 15,16 However, because latent classes were identified using an unbiased, data-driven approach, the results are comparable to those of previous reports that use data from a RCT, allowing potential causal inferences regarding treatment effects and potential similarities to previously reported subphenotypes. 31

Original manuscripts

Calfee CS, Delucchi KL, Sinha P, Matthay MA, Hackett J, Shankar-Hari M, et al. Acute respiratory distress syndrome subphenotypes and differential response to simvastatin: secondary analysis of a randomised controlled trial. Lancet Respir Med 2018;6:691–8.

Santhakumaran S, Gordon A, Prevost AT, O’Kane C, McAuley DF, Shankar-Hari M. Heterogeneity of treatment effect by baseline risk of mortality in critically ill patients: re-analysis of three recent sepsis and ARDS randomised controlled trials. Crit Care 2019;23:156.

Reviews and editorials

Shankar-Hari M, McAuley DF. Acute respiratory distress syndrome phenotypes and identifying treatable traits. The dawn of personalized medicine for ARDS. Am J Respir Crit Care Med 2017;195:280–1.

Shankar-Hari M, McAuley DF. Divide and conquer: identifying acute respiratory distress syndrome subphenotypes. Thorax 2017;72:867–9.

Shankar-Hari M, Fan E, Ferguson ND. Acute respiratory distress syndrome (ARDS) phenotyping. Intensive Care Med 2019;45:516–19.

Poster presentation

Ferris P, Boyle A, Conlon J, Gordon AC, Shankar-Hari M, O’Kane C, McAuley D. Baseline NT-proBNP Predicts outcOme and Treatment Response to Statin Therapy in Patients with ARDS. American Thoracic Society 2019 International Conference, Dallas, TX, USA, 17–22 May 2019.