LABOR Study

Author

Russell Lewis

Published

September 2, 2024

Background

This is the webpage summarizes ongoing analysis for the LABOR study. The figures are preliminary and in many cases I have not changed the variable indicators/legends from the original dataset (I don’t have a code book of what some the numbers represent). The figures can be aesthetically improved and reformatted before final publication. This analysis is just an initial attempt to understand the dataset and confirm a data analysis approach for the final manuscript. The working hypothesis of the analysis is that L-AMB is non-inferior to voriconazole/isavuconazole (triazole-based therapy) for the treatment of invasive aspergillosis using using Day + 90 all-cause mortality as the primary endpoint.

To save time, I did not recreate your “Table 1” analysis you forwarded to me, but used it in consideration of preparing the analysis.

Primary endpoint

The first few graphs are an attempt to just visualize the data. I started by plotting the Kaplan-Meier Curves for the entire dataset stratified by whether patients received triazoles or L-AMB.

Code

data<-read_xlsx("LABOR_definitivo.xlsx")
fit <- survfit(Surv(IA_to_90death_COX, death_3m) ~ terapia_antifungina_categorica , data = data)

Code

Kmgraph <-ggsurvplot(
  fit,
  data = data,
  size = 1,                 # change line size
  palette =
    c("red", "#2E9FDF"),# custom color palettes
  conf.int = TRUE,          # Add confidence interval
  pval = TRUE,              # Add p-value
  risk.table = TRUE,        # Add risk table
  risk.table.col = "strata",# Risk table color by groups

  legend.labs =
    c("Triazole", "Liposomal AMB"),    # Change legend labels
  risk.table.height = 0.25, # Useful to change when you have multiple groups
  ggtheme = theme_bw()      # Change ggplot2 theme
)

Kmgraph

Figure 1: Kaplan-Meier analysis of 90-day survival rates

These data demonstrate that overall, there is no significant differences (using Log-Rank test) in the rate of 90-day mortality between patients who received L-AMB versus triazoles, although there seems to be a trend favoring the triazoles.

To confirm that the rates of mortality are similarly distributed between the two groups I created population time plots using the casebase package¹. each grey line represents the time course of a single patients (the graph is arranged from shortest to longest) with red dots indicating the time of mortality. The distribution in time to death is relatively similar between the two groups although the highest incidence density of mortality appears in the first 25 days for both groups. This is consistent with using Day +42 mortality as a more indicative endpoint for fungal infection-specific mortality

Code

library (data.table)

pt_fungal <- popTime(data = data, time="IA_to_90death_COX", event="death_3m")
class(pt_fungal)

[1] "popTime"    "data.table" "data.frame"

Code

plot(popTime(pt_fungal, exposure = "terapia_antifungina_categorica"))

Figure 2: Population time-plot of 90-day mortality, 1= Triazole; 2=Liposomal AMB; Each gray line represent 1 subject time in study. Red dots indicate time of death

I noticed in your Table 1 there was significant differences in the switching to other antifungals between L-AMB and voriconazole patients. I reset the data to day +90 for the three key time variables (LOS in hospital, time to switch and time to death by 90 days) and compared the cumulative incidence curves. We can add length of ICU stay (LIS) later if desired although it will make the graph harder to read. Overall the mortality CID curves and hospital discharge curves are similar. However These data suggest that more than 2/3 of patients who received L-AMB were switched to an alternative antifungal within two weeks. This is consistent with empirical use and best clinical practice for use of the drug to avoid nephrotoxicity,² but creates a bias in trying to demonstrate “equivalence” of the two treatment approaches. Therefore we will mostly be comparing a strategy of starting two weeks with liposomal amphotericin B.

Analysis option

Consider a table that describes patients (n=38) with aspergillosis who did not switch from liposomal amphotericin B

Code

library (readxl)
library (ggplot2)
cif_LAMB <- read_excel("cif_LAMB.xlsx")
cif_voriconazole <- read_excel("cif_voriconazole.xlsx")

graph1<-ggplot () +
  geom_area(data = cif_LAMB, aes(x = time, y = Surv), fill = "red", alpha = 0.5) +
  geom_area(data = cif_LAMB, aes(x = time, y = Dschg), fill = "green", alpha = 0.5) +
  geom_area(data = cif_LAMB, aes(x = time, y = Switch), fill = "blue", alpha = 0.5) +
  labs(title = "Amphotericin B", x = "Days since diagnosis", y = "Incidence") + theme_minimal()

graph2<-ggplot () +
  geom_area(data = cif_voriconazole, aes(x = time, y = Surv), fill = "red", alpha = 0.5) +
  geom_area(data = cif_voriconazole, aes(x = time, y = Dschg), fill = "green", alpha = 0.5) +
  geom_area(data = cif_voriconazole, aes(x = time, y = Switch), fill = "blue", alpha = 0.5) +
  labs(title = "Voriconazole", x = "Days since diagnosis", y = "Incidence") + theme_minimal()

library("cowplot")
plot_grid(graph1, graph2,
          ncol = 2, nrow = 1)

Figure 3: Cumulative incidence function of death (red), discharge (green) and medication switch (blue) between the two study cohorts

Indeed, if we go back and compare the Kaplan-Meier curves for liposomal AMB and voriconazole in patients who did not switch therapies, we see a significant difference in mortality by the KM curves. However, I recognize this analysis may be biased against L-AMB.

Code

Kmgraph2 <-ggsurvplot(
  fit2,
  data = filtered_data,
  size = 1,                 
  palette =
    c("red", "#2E9FDF"),
  conf.int = TRUE,         
  pval = TRUE,           
  risk.table = TRUE,      
  risk.table.col = "strata",
  legend.labs =
    c("Triazole", "Liposomal AMB"),    # Change legend labels
  risk.table.height = 0.25, # Useful to change when you have multiple groups
  ggtheme = theme_bw()      # Change ggplot2 theme
)

Kmgraph2

Figure 4: Kaplan-Meier analysis in patients who did not switch therapy

A comparison of the specific switch in therapy shows marked differences in the types of molecules that were substituted (I do not know specifically which drugs these numerical codes represent).

Code

vori_change <-ggplot(vori_df, aes(x = tp_modif_molecola, y = percent, fill = tp_modif_molecola)) +
  geom_bar(stat = "identity") +
  labs(title = "Triazole", x = "Molecule", y = "Percent") +
  theme_minimal() + scale_y_continuous(labels = scales::percent_format(scale = 1), limits = c(0, 100)) +
   guides(fill = "none")

lamb_change <-ggplot(lamb_df, aes(x = tp_modif_molecola, y = percent, fill = tp_modif_molecola)) +
  geom_bar(stat = "identity") +
  labs(title = "Liposomal AMB", x = "Molecule", y = "Percent") +
  theme_minimal() + scale_y_continuous(labels = scales::percent_format(scale = 1), limits = c(0, 100)) +
   guides(fill = "none")

library("cowplot")
plot_grid(vori_change, lamb_change,
          ncol = 2, nrow = 1)

Figure 5: Therapy changes by type of antifungal molecule

The reasons were switching were also different, especially reason #1 (nephrotoxicity or step-down to oral therapy?). I noticed that there is no variable in the data set for nephrotoxicity?

Code

#| label: fig-change 2
#| fig-cap: "Therapy changes by reason for change"
#| warning: false
vori_change2 <-ggplot(vori_df2, aes(x = tp_modif_motivo, y = percent, fill = tp_modif_motivo)) +
  geom_bar(stat = "identity") +
  labs(title = "Triazole", x = "Molecule", y = "Percent") +
  theme_minimal() + scale_y_continuous(labels = scales::percent_format(scale = 1), limits = c(0, 100)) +
   guides(fill = "none")

lamb_change2 <-ggplot(lamb_df2, aes(x = tp_modif_motivo, y = percent, fill = tp_modif_motivo)) +
  geom_bar(stat = "identity") +
  labs(title = "Liposomal AMB", x = "Molecule", y = "Percent") +
  theme_minimal() + scale_y_continuous(labels = scales::percent_format(scale = 1), limits = c(0, 100)) +  
   guides(fill = "none")

library("cowplot")
plot_grid(vori_change2, lamb_change2,
          ncol = 2, nrow = 1)

Cox regression: Survival

I then examined possible multivariable models that could potentially allow for comparisons accounting for differences in patient risk factores for 90-day mortality. Examining your univariate analysis, I included several key risk factors that would be expected or were shown in your analysis to be associated with a higher risk of 90-mortality. For the moment I did not include biomarkers (I will come back to these later). Only three variables emerged as independent risk factors for 90-mortality: CCI, fr_VAPA, and terapia_antifungina_categorica. Comparison of residuals versus linear predictors did not identify a need for modeling these as time-varying covariates (also shown by the chi-square statistics vs. residuals).

Code

res.cox <- coxph(Surv(IA_to_90death_COX, death_3m) ~ eta + CCI + fr_SOT + fr_HSCT + fr_malatt_ematol +  fr_neutropenia + fr_steroide_cronico + fr_VAPA + fr_degenzaTI_sup14gg + profilassi_antifungina + terapia_antifungina_categorica, data = data_tibble)
summary(res.cox)

Call:
coxph(formula = Surv(IA_to_90death_COX, death_3m) ~ eta + CCI + 
    fr_SOT + fr_HSCT + fr_malatt_ematol + fr_neutropenia + fr_steroide_cronico + 
    fr_VAPA + fr_degenzaTI_sup14gg + profilassi_antifungina + 
    terapia_antifungina_categorica, data = data_tibble)

  n= 401, number of events= 171 

                                    coef exp(coef)  se(coef)      z Pr(>|z|)   
eta                             0.008539  1.008576  0.007938  1.076  0.28203   
CCI                             0.097698  1.102629  0.037363  2.615  0.00893 **
fr_SOT                         -0.471556  0.624031  0.291787 -1.616  0.10607   
fr_HSCT                        -0.308537  0.734521  0.296504 -1.041  0.29807   
fr_malatt_ematol               -0.084158  0.919286  0.208580 -0.403  0.68660   
fr_neutropenia                 -0.115370  0.891037  0.264951 -0.435  0.66324   
fr_steroide_cronico            -0.093865  0.910405  0.227392 -0.413  0.67976   
fr_VAPA                         0.337957  1.402080  0.176005  1.920  0.05484 . 
fr_degenzaTI_sup14gg            0.164585  1.178904  0.254167  0.648  0.51728   
profilassi_antifungina          0.150528  1.162448  0.214695  0.701  0.48323   
terapia_antifungina_categorica  0.535144  1.707695  0.192275  2.783  0.00538 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

                               exp(coef) exp(-coef) lower .95 upper .95
eta                               1.0086     0.9915    0.9930     1.024
CCI                               1.1026     0.9069    1.0248     1.186
fr_SOT                            0.6240     1.6025    0.3522     1.106
fr_HSCT                           0.7345     1.3614    0.4108     1.313
fr_malatt_ematol                  0.9193     1.0878    0.6108     1.384
fr_neutropenia                    0.8910     1.1223    0.5301     1.498
fr_steroide_cronico               0.9104     1.0984    0.5830     1.422
fr_VAPA                           1.4021     0.7132    0.9930     1.980
fr_degenzaTI_sup14gg              1.1789     0.8482    0.7164     1.940
profilassi_antifungina            1.1624     0.8603    0.7632     1.771
terapia_antifungina_categorica    1.7077     0.5856    1.1715     2.489

Concordance= 0.631  (se = 0.021 )
Likelihood ratio test= 39.26  on 11 df,   p=5e-05
Wald test            = 38.77  on 11 df,   p=6e-05
Score (logrank) test = 39.04  on 11 df,   p=5e-05

Code

## Test of proportaional hazards assumption
#| label: fig-cox2
#| fig-cap: "Tests of proportional hazards assumption "
#| warning: false
test.ph <- cox.zph(res.cox)
test.ph

                                  chisq df    p
eta                            0.155110  1 0.69
CCI                            0.317143  1 0.57
fr_SOT                         0.000987  1 0.97
fr_HSCT                        0.388241  1 0.53
fr_malatt_ematol               0.687075  1 0.41
fr_neutropenia                 0.349761  1 0.55
fr_steroide_cronico            1.446334  1 0.23
fr_VAPA                        0.054621  1 0.82
fr_degenzaTI_sup14gg           0.661913  1 0.42
profilassi_antifungina         0.121841  1 0.73
terapia_antifungina_categorica 0.321447  1 0.57
GLOBAL                         5.681352 11 0.89

Code

ggcoxdiagnostics(res.cox, type = "dfbeta" , linear.predictions = TRUE)

Figure 6: Multivariable Cox regression analysis of risk factors for 90-day mortality

Code

ggforest(
 res.cox,
  data = data,
  main = "Hazard ratio",
  cpositions = c(0.02, 0.22, 0.4),
  fontsize = 0.7,
  refLabel = "reference",
  noDigits = 2
)

Figure 7: Plot of Cox regression hazard ratios

When survival curves for triazole versus L-AMB therapy were adjusted for the two independent predictors of 90-day mortality (CCI and VAPA) we see that higher survival rates were predicted for patients in the triazole (group 1) versus liposomal AMB (group 2) treatment. The lower graph shows the estimated mortality difference over time (between 5-10%) difference with the 95% CI climbing above zero by approximately day 40.

Code

library (survival)
data$terapia_antifungina_categorica <- as.factor(data$terapia_antifungina_categorica)
outcome_model <- survival::coxph(Surv(IA_to_90death_COX, death_3m) ~ CCI + fr_VAPA + terapia_antifungina_categorica, data = data, x=TRUE)
library (adjustedCurves)
adjsurv <- adjustedsurv(data=data,
                        variable="terapia_antifungina_categorica",
                        ev_time="IA_to_90death_COX",
                        event="death_3m",
                        method="direct",
                        outcome_model=outcome_model,
                        conf_int=TRUE)
head(adjsurv$adj)

  time      surv group          se  ci_lower  ci_upper
1    0 1.0000000     1 0.000000000 1.0000000 1.0000000
2    2 0.9954973     1 0.003164160 0.9892957 1.0000000
3    3 0.9796893     1 0.006483467 0.9669819 0.9923966
4    4 0.9728777     1 0.007457499 0.9582613 0.9874941
5    5 0.9592456     1 0.009241315 0.9411329 0.9773582
6    6 0.9501513     1 0.010256284 0.9300493 0.9702532

Code

adjusted_curve_diff(adjsurv, times=14, conf_int=TRUE)

  time       diff         se     ci_lower  ci_upper    p_value
1   14 0.06026799 0.03580852 -0.009915409 0.1304514 0.09236288

Code

adjusted_curve_diff(adjsurv, times=30, conf_int=TRUE)

  time       diff         se     ci_lower  ci_upper    p_value
1   30 0.08563529 0.04457766 -0.001735308 0.1730059 0.05472714

Code

adjusted_curve_diff(adjsurv, times=42, conf_int=TRUE)

  time      diff         se    ci_lower  ci_upper    p_value
1   42 0.0983702 0.04864434 0.003029043 0.1937114 0.04315227

Code

adjusted_curve_diff(adjsurv, times=60, conf_int=TRUE)

  time      diff         se    ci_lower  ci_upper  p_value
1   60 0.1116029 0.05194426 0.009794052 0.2134118 0.031673

Code

adjusted_curve_diff(adjsurv, times=90, conf_int=TRUE)

  time      diff         se   ci_lower  ci_upper    p_value
1   90 0.1186066 0.05351633 0.01371651 0.2234967 0.02667308

Code

adjusted_rmst(adjsurv, to=42, conf_int=TRUE)

  to group     rmst
1 42     1 34.92457
2 42     2 32.27469

Code

library (ggthemes)
library (ggplot2)
coxplot1 <-plot(adjsurv, conf_int=TRUE, custom_colors=c("red", "blue"), gg_theme=theme_bw(), risk_table_theme=theme_classic(),
     legend.position="top", xlab="Time (days)") +
  ylim(0, 1)
coxplot2 <-plot_curve_diff(adjsurv, conf_int=TRUE, color="blue")

library("cowplot")
plot_grid(coxplot1, coxplot2,
          ncol = 2, nrow = 1)

Figure 8: Kaplan-Meier survival curves adjusted for 90-day mortality risk factors. Top graph: Group 1= triazoles, Group 2= L-AMB; Lower graph- Estimated survival differences over time between triazole vs. L-AMB treated patients.

Biomarkers

Galactomannan

I examined the relationship between serum or BAL galactomannan index versus 90-day mortality by logistic regression. A weak relationship between increasing serum galactomannan and survival was evident (left panel); but no relationship was evident by BAL GM index and mortality (right panel).

Code

library (jtools)
library (cowplot)

model1 <- glm(death_3m ~ GM_siero_valore, data = data, family = binomial)
summ(model1)

MODEL INFO:
Observations: 87 (314 missing obs. deleted)
Dependent Variable: death_3m
Type: Generalized linear model
  Family: binomial 
  Link function: logit 

MODEL FIT:
χ²(1) = 3.14, p = 0.08
Pseudo-R² (Cragg-Uhler) = 0.05
Pseudo-R² (McFadden) = 0.03
AIC = 121.46, BIC = 126.39 

Standard errors:MLE
----------------------------------------------------
                         Est.   S.E.   z val.      p
--------------------- ------- ------ -------- ------
(Intercept)             -0.48   0.34    -1.41   0.16
GM_siero_valore          0.22   0.13     1.70   0.09
----------------------------------------------------

Code

data$predicted_mortality <- predict(model1, newdata = data, type = "response")
nrow(data)

[1] 401

Code

length(predict(model1, newdata = data, type = "response"))

[1] 401

Code

effect1 <-effect_plot(model1, pred = GM_siero_valore, interval = TRUE, plot.points = TRUE, 
            jitter = 0.05)

model2 <- glm(death_3m ~ GM_BAL_valore, data = data, family = binomial)
summary (model2)


Call:
glm(formula = death_3m ~ GM_BAL_valore, family = binomial, data = data)

Coefficients:
              Estimate Std. Error z value Pr(>|z|)
(Intercept)   -0.12940    0.24355  -0.531    0.595
GM_BAL_valore -0.02891    0.05935  -0.487    0.626

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 457.06  on 332  degrees of freedom
Residual deviance: 456.82  on 331  degrees of freedom
  (68 observations deleted due to missingness)
AIC: 460.82

Number of Fisher Scoring iterations: 3

Code

data$predicted_mortality2 <- predict(model2, newdata = data, type = "response")

nrow(data)

[1] 401

Code

length(predict(model2, newdata = data, type = "response"))

[1] 401

Code

effect2 <-effect_plot(model2, pred = GM_BAL_valore, interval = TRUE, plot.points = TRUE, 
            jitter = 0.05)

library("cowplot")
plot_grid(effect1, effect2,
          ncol = 2, nrow = 1)

Figure 9: Relationship of galactomannan index in serum or BAL versus 90-day mortality

PCR

I also compared log10 transformed PCR copy data and 90-day mortality. No clear relationship between PCR copy number and mortality was evident.

Code

#| label: fig-gm
#| fig-cap: "Relationship of log-transformed serum PCR copies versus 90-day mortality"
#| warning: false
library (jtools)
library (cowplot)



data$log10_variable <- log10(data$copie_PCR)

model3 <- glm(death_3m ~log10_variable , data = data, family = binomial)
summ(model3)

MODEL INFO:
Observations: 43 (358 missing obs. deleted)
Dependent Variable: death_3m
Type: Generalized linear model
  Family: binomial 
  Link function: logit 

MODEL FIT:
χ²(1) = 0.08, p = 0.78
Pseudo-R² (Cragg-Uhler) = 0.00
Pseudo-R² (McFadden) = 0.00
AIC = 61.64, BIC = 65.16 

Standard errors:MLE
---------------------------------------------------
                        Est.   S.E.   z val.      p
-------------------- ------- ------ -------- ------
(Intercept)            -0.87   1.64    -0.53   0.59
log10_variable          0.13   0.46     0.28   0.78
---------------------------------------------------

Code

data$predicted_mortality <- predict(model3, newdata = data, type = "response")
nrow(data)

[1] 401

Code

length(predict(model3, newdata = data, type = "response"))

[1] 401

Code

effect_plot(model3, pred =log10_variable, interval = TRUE, plot.points = TRUE, 
                      jitter = 0.05)

Summary

In this retrospective analysis, patients who initially received triazole-based regimens are predicted to have improved survival versus liposomal AMB
- This difference was more evident after controlling for CCI and VAPA
- Comparisons are complicated by frequent switching of therapy, therefore differences at time of switch and earlier timepoints could also be investigated or reported (e.g., day 30, day 42)
- Stratified analysis based on switch may be informative but will reduce power of analysis
There was a modest relationship between serum GM index and 90-day mortality, but no clear relationship between BAL GM index or serum PCR copies and mortality
Overall survival curves were consistent with those previously reported in the literature

References

Bhatnagar, S., Turgeon, M., Islam, J., Saarela, O. & Hanley, J. Casebase: An alternative framework for survival analysis and comparison of event rates. 14, (2022).

Stanzani, M. et al. Retrospective cohort analysis of liposomal amphotericin b nephrotoxicity in patients with hematological malignancies. Antimicrobial Agents and Chemotherapy 61, (2017).