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)summary(survfit(Surv(IA_to_90death_COX, death_3m) ~ terapia_antifungina_categorica, data = data), times =42)
Call: survfit(formula = Surv(IA_to_90death_COX, death_3m) ~ terapia_antifungina_categorica,
data = data)
terapia_antifungina_categorica=1
time n.risk n.event survival std.err lower 95% CI
42.0000 208.0000 88.0000 0.7027 0.0266 0.6525
upper 95% CI
0.7567
terapia_antifungina_categorica=2
time n.risk n.event survival std.err lower 95% CI
42.0000 69.0000 36.0000 0.6571 0.0463 0.5723
upper 95% CI
0.7545
Code
summary(survfit(Surv(IA_to_90death_COX, death_3m) ~ terapia_antifungina_categorica, data = data), times =90)
Call: survfit(formula = Surv(IA_to_90death_COX, death_3m) ~ terapia_antifungina_categorica,
data = data)
terapia_antifungina_categorica=1
time n.risk n.event survival std.err lower 95% CI
90.0000 177.0000 122.0000 0.5878 0.0286 0.5344
upper 95% CI
0.6467
terapia_antifungina_categorica=2
time n.risk n.event survival std.err lower 95% CI
90.0000 58.0000 49.0000 0.5333 0.0487 0.4460
upper 95% CI
0.6378
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
Changes in 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).
Figure 7: Therapy changes by type of antifungal molecule
Matteo: I noticed there were actually 8 reasons in the original dataset as motives for changing therapy, but in your email 03.12.2024 only 7 reasons were given. What is the 8th reason?
I then examined possible multivariable models that could potentially allow for comparisons accounting for differences in patient risk factors 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 than 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 category of antifungal therapy remain in the model. However, marked changes in the HR for fr_degenzaTI_sup14gg and ricovero_TI suggest a problem with colinearity 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
“Dove_diagnosi=3” was suggested as the best indicator value for ICU status. So I created a new dummay variable for icu “new_icu” and re-ran the cox-regression. ICU status, CCI and the type of antifungal therapy remained as independent predict variables of 90-day mortality.
Code
library(dplyr)data <- data %>%mutate(new_icu =ifelse(dove_diagnosi ==3, 1, 0))res.cox3 <-coxph(Surv(IA_to_90death_COX, death_3m) ~ eta + CCI + fr_SOT + fr_HSCT + fr_malatt_ematol + fr_neutropenia + fr_steroide_cronico + fr_VAPA + new_icu + profilassi_antifungina + terapia_antifungina_categorica, data = data)summary(res.cox3)
Very few patients with VAPA received L-AMB, therefore it is not surprising 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 14: 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 15: Comparison of Kaplan-Meier Survival curves in patients with VAPA
Predicted Cox regression survival curves by treatment type
When survival curves for triazole versus L-AMB therapy were adjusted for the two independent predictors of 90-day mortality (CCI and ICU_new) 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 ( nearly 10%) difference with the 95% CI climbing suprassing above zero within the first 25 days,
Code
library (survival)library(dplyr)data_icu <- data %>%mutate(new_icu =ifelse(dove_diagnosi ==3, 1, 0))data_icu$terapia_antifungina_categorica <-as.factor(data_icu$terapia_antifungina_categorica)outcome_model <- survival::coxph(Surv(IA_to_90death_COX, death_3m) ~ CCI + new_icu + terapia_antifungina_categorica, data = data_icu, x=TRUE)library (adjustedCurves)adjsurv <-adjustedsurv(data=data_icu,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 16: 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.
Cox regression predicted survival
When survival curves for triazole versus L-AMB therapy are analyzed in relation to icu_new and predictors for mortalty (icu_new, CCI) we 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)library(dplyr)data_icu <- data %>%mutate(new_icu =ifelse(dove_diagnosi ==3, 1, 0))data_icu$terapia_antifungina_categorica <-as.factor(data_icu$terapia_antifungina_categorica)outcome_model2 <- survival::coxph(Surv(IA_to_90death_COX, death_3m) ~ CCI + new_icu + fr_SOT + terapia_antifungina_categorica, data = data_icu, x=TRUE)library (adjustedCurves)adjsurv2 <-adjustedsurv(data=data_icu,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 17: 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 18: 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)
Voriconazole TDM and Survival
First sample
In the analysis of the first sample TDM and outcome, no significant relationship was observed by linear regression between the first voriconazole TDM result and 90-day mortality.
Code
filter_tdm <- data %>%filter (prima_tdm_vorico >0.05)ggplot(filter_tdm, aes(x = prima_tdm_vorico)) +geom_histogram(binwidth =2, fill ="skyblue", color ="black", alpha =0.7) +labs(title ="First voriconazole TDM", x ="Concentrations (mg/L)", y ="Frequency") +theme_minimal()
Code
filter_tdm <- data %>%filter (prima_tdm_vorico >0.05)model4 <-glm(death_3m ~ prima_tdm_vorico, data = filter_tdm, family = binomial)summ(model4)
MODEL INFO:
Observations: 191
Dependent Variable: death_3m
Type: Generalized linear model
Family: binomial
Link function: logit
MODEL FIT:
χ²(1) = 0.46, p = 0.50
Pseudo-R² (Cragg-Uhler) = 0.00
Pseudo-R² (McFadden) = 0.00
AIC = 261.11, BIC = 267.61
Standard errors:MLE
-----------------------------------------------------
Est. S.E. z val. p
---------------------- ------- ------ -------- ------
(Intercept) -0.51 0.23 -2.25 0.02
prima_tdm_vorico 0.03 0.04 0.67 0.50
-----------------------------------------------------
Code
filter_tdm$predicted_mortality <-predict(model4, newdata = filter_tdm, type ="response")nrow(filter_tdm)
[1] 191
Code
length(predict(model4, newdata = filter_tdm, type ="response"))
Relationship between second voriconazole TDM and outcome
In the analysis of the second voriconazole TDM and outcome, a significant relationship was observed by logistic regression for TDM concentrations and 90-survival. Hower it appears that patients with very high voriconazole TDM concentrations had higher mortality, which may be a reflection of organ failure.
Code
filter_tdm2 <- data %>%filter (seconda_tdm_vorico >0.05)ggplot(filter_tdm2, aes(x = seconda_tdm_vorico)) +geom_histogram(binwidth =2, fill ="skyblue", color ="black", alpha =0.7) +labs(title ="Second voriconazole TDM", x ="Concentrations (mg/L)", y ="Frequency") +theme_minimal()
Code
filter_tdm2 <- data %>%filter (seconda_tdm_vorico >0.05)model5 <-glm(death_3m ~ seconda_tdm_vorico, data = filter_tdm2, family = binomial)summ(model5)
MODEL INFO:
Observations: 147
Dependent Variable: death_3m
Type: Generalized linear model
Family: binomial
Link function: logit
MODEL FIT:
χ²(1) = 7.54, p = 0.01
Pseudo-R² (Cragg-Uhler) = 0.07
Pseudo-R² (McFadden) = 0.04
AIC = 188.65, BIC = 194.64
Standard errors:MLE
-------------------------------------------------------
Est. S.E. z val. p
------------------------ ------- ------ -------- ------
(Intercept) -1.26 0.32 -3.98 0.00
seconda_tdm_vorico 0.15 0.06 2.66 0.01
-------------------------------------------------------
Code
filter_tdm2$predicted_mortality <-predict(model5, newdata = filter_tdm2, type ="response")nrow(filter_tdm2)
[1] 147
Code
length(predict(model4, newdata = filter_tdm2, type ="response"))
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 for several confounding factors
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 interest 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
Changelog
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
11.12.2024: Request to use “dove_diagnosi_3” for ICU variable in the Cox regression -Completed
11.12.2024: I will also examine relationship of TDM endpoints with outcomes - in progress
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).
