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 sizepalette =c("red", "#2E9FDF"),# custom color palettesconf.int =TRUE, # Add confidence intervalpval =TRUE, # Add p-valuerisk.table =TRUE, # Add risk tablerisk.table.col ="strata",# Risk table color by groupslegend.labs =c("Triazole", "Liposomal AMB"), # Change legend labelsrisk.table.height =0.25, # Useful to change when you have multiple groupsggtheme =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.
Comparison of ICU versus non- ICU patients
Interestingly, differences between l_ABM and triazoles is not observed in patients who were in the ICU, but patients receiving triazoles who were not in the ICU had lower mortality rates
Code
# Filter data for patients in or not in ICUICU_data <-filter(data, ricovero_TI =="2")nonicu_data <-filter(data, ricovero_TI =="1")
Code
fiticu <-survfit(Surv(IA_to_90death_COX, death_3m) ~ terapia_antifungina_categorica , data = ICU_data)Kmgraph_ICU <-ggsurvplot( fiticu,data = ICU_data,size =1, # change line sizepalette =c("red", "#2E9FDF"),# custom color palettesconf.int =TRUE, # Add confidence intervalpval =TRUE, # Add p-valuerisk.table =TRUE, # Add risk tablerisk.table.col ="strata",# Risk table color by groupslegend.labs =c("Triazole", "Liposomal AMB"), # Change legend labelsrisk.table.height =0.25, # Useful to change when you have multiple groupsggtheme =theme_bw() # Change ggplot2 theme)Kmgraph_ICU
Figure 2: Kaplan-Meier analysis of 90-day survival rates in ICU patients
Code
fitnonicu <-survfit(Surv(IA_to_90death_COX, death_3m) ~ terapia_antifungina_categorica , data = nonicu_data)Kmgraph_nonICU <-ggsurvplot( fitnonicu,data = nonicu_data,size =1, # change line sizepalette =c("red", "#2E9FDF"),# custom color palettesconf.int =TRUE, # Add confidence intervalpval =TRUE, # Add p-valuerisk.table =TRUE, # Add risk tablerisk.table.col ="strata",# Risk table color by groupslegend.labs =c("Triazole", "Liposomal AMB"), # Change legend labelsrisk.table.height =0.25, # Useful to change when you have multiple groupsggtheme =theme_bw() # Change ggplot2 theme)Kmgraph_nonICU
Figure 3: Kaplan-Meier analysis of 90-day survival rates in non-ICU patients
Mortality rate distribution
To confirm that the rates of mortality are similarly distributed between the two groups I created population time plots using the casebase package1. 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
Figure 4: 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,2 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 5: Cumulative incidence function of survival (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 labelsrisk.table.height =0.25, # Useful to change when you have multiple groupsggtheme =theme_bw() # Change ggplot2 theme)Kmgraph2
Figure 6: 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 7: 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: falsevori_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). Note this analysis included ICU stay longer thant 14 days.
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)
Based on our discussions I re-added any ICU stay in addition to 14 days stay in the ICU in the model. The results are somewhat different as VAPA drops from the model, SOT, ICU stay, CCI and cateogry of antifungal therapy remain in the model
## Test of proportional hazards assumption#| label: fig-cox3#| fig-cap: "Tests of proportional hazards assumption with revised model "#| warning: falsetest.ph2 <-cox.zph(res.cox2)test.ph2
Very few patients with VAPA received L-AMB, therefore it is not surpsiging that no differences are seen when patients are stratified by VAPA
Code
novapa_data <-filter(data_tibble, fr_VAPA =="0")vapa_data <-filter(data_tibble, fr_VAPA =="1")fit_novapa <-survfit(Surv(IA_to_90death_COX, death_3m) ~ terapia_antifungina_categorica , data = novapa_data)fit_vapa <-survfit(Surv(IA_to_90death_COX, death_3m) ~ terapia_antifungina_categorica , data = vapa_data)Kmgraph_novapa <-ggsurvplot( fit_novapa,data = novapa_data,size =1, # change line sizepalette =c("red", "#2E9FDF"),# custom color palettesconf.int =TRUE, # Add confidence intervalpval =TRUE, # Add p-valuerisk.table =TRUE, # Add risk tablerisk.table.col ="strata",# Risk table color by groupslegend.labs =c("Triazole", "Liposomal AMB"), # Change legend labelsrisk.table.height =0.25, # Useful to change when you have multiple groupsggtheme =theme_bw() # Change ggplot2 theme)Kmgraph_novapa
Figure 12: Comparison of Kapla-Meier Survival curves in patients without VAPA
Code
Kmgraph_vapa <-ggsurvplot( fit_vapa,data = vapa_data,size =1, # change line sizepalette =c("red", "#2E9FDF"),# custom color palettesconf.int =TRUE, # Add confidence intervalpval =TRUE, # Add p-valuerisk.table =TRUE, # Add risk tablerisk.table.col ="strata",# Risk table color by groupslegend.labs =c("Triazole", "Liposomal AMB"), # Change legend labelsrisk.table.height =0.25, # Useful to change when you have multiple groupsggtheme =theme_bw() # Change ggplot2 theme)Kmgraph_vapa
Figure 13: Comparison of Kaplan-Meier Survival curves in patients with VAPA
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)
Figure 14: 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.
When survival curves for triazole versus L-AMB therapy are repeated to including the additional predictors for mortalty (ICU, CCI, and SOT) we still see 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_model2 <- survival::coxph(Surv(IA_to_90death_COX, death_3m) ~ CCI + ricovero_TI + fr_SOT + terapia_antifungina_categorica, data = data, x=TRUE)library (adjustedCurves)adjsurv2 <-adjustedsurv(data=data,variable="terapia_antifungina_categorica",ev_time="IA_to_90death_COX",event="death_3m",method="direct",outcome_model=outcome_model2,conf_int=TRUE)head(adjsurv2$adj)
Figure 15: Kaplan-Meier survival curves adjusted for 90-day mortality risk factors including additional 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"))
Figure 16: 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: falselibrary (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 is still evident after controlling
Comparisons are complicated by frequent switching of therapy- i.e. patients who are clinically resposing tp L-AMB are probably switched to triazoles
Stratified analysis is of iterest but reduced the power; so no differences are seen between regimens
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
Follow-up analysis after September 6th Teams call
Separation of ICU vs non-ICU -Completed
Inclusion of ICU in MV model -Completed
Separation of VAPA -Completed
Multistate analysis by type of therapy? - not attempted because of power issues
References
1.
Bhatnagar, S., Turgeon, M., Islam, J., Saarela, O. & Hanley, J. Casebase: An alternative framework for survival analysis and comparison of event rates. 14, (2022).
