Upload New File

ICP_burden_finalfigs.R

+# Creates figures for the manuscript
+
+
+rm(list=ls())
+
+library(ggplot2)
+library(gridExtra)
+library(dplyr)
+library(tidyr)
+library(heatwaveR)
+library(ggpubr)
+library(ggsci)
+library(RColorBrewer)
+library(grid)
+
+nas = c(NA, "")
+f <- "C:/Users/Cecilia KI/Dropbox/Cecilia/ICP burden/ICP_tempdata/200625/"
+jet.colors <- colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F", "yellow", "#FF7F00", "red"))#, "#7F0000"))
+load(file=paste0(f,"PTD.RData"))
+load(paste0(f,"Corr.list.RData"))
+Corr_bootstraps.GOSE <- read.csv(paste0(f,"Corr_bootstraps.GOSE6mo.df.10.csv"), stringsAsFactors = FALSE)
+Autoregtimes <- read.csv(paste0(f,"Autoregtimes.csv"), stringsAsFactors = FALSE, na.strings = nas)
+
+ms <- read.csv("C:/Users/Cecilia KI/Dropbox/Cecilia/ICP burden/ms.new_200625.csv", na.strings = nas, stringsAsFactors = FALSE)
+
+ms <- ms %>% mutate(datetime = as.POSIXct(datetime),
+                    gupi = as.factor(gupi),
+                    GOS = as.factor(GOS),
+                    GOSE6mo = as.factor(GOSE6mo))
+
+ms <- ms %>% mutate(icp= ifelse(is.na(icp), icp.evd, icp))
+
+ms <- ms %>% group_by(gupi) %>% mutate(monitor.time = as.numeric(difftime(max(datetime), min(datetime), units = "days")))
+
+Baseline <- read.csv("C:/Users/Cecilia KI/Dropbox/Cecilia/ICP burden/Baseline.new_200625.csv", stringsAsFactors = FALSE, na.strings = nas) #Created in script ICP_burden_hypothermia_preparation_190607.csv
+
+Baseline <- Baseline %>% mutate(Mort90d = as.factor(ifelse(as.Date(DeathDate) <= as.Date("1970-01-01") + 90, 1,0)),
+                                Mort6mo = as.factor(ifelse(as.Date(DeathDate) <= as.Date("1970-01-01") + 180, 1, 0)),
+                                Outcome = as.factor(ifelse(as.integer(GOSE6mo) <= 4, "UNFAV",
+                                                           ifelse(as.integer(GOSE6mo) > 4, "FAV", NA))),
+                                GOSE6mo = as.factor(GOSE6mo))
+Baseline$Mort90d[is.na(Baseline$Mort90d)] <- 0
+Baseline$Mort6mo[is.na(Baseline$Mort6mo)] <- 0
+Baseline <- Baseline %>% mutate(Mort6mo = ifelse(Mort6mo == 1, "Yes", "No"))
+
+events_all <- read.csv(paste0(f,"events_all_200625.csv"), stringsAsFactors = FALSE, na.strings = nas)
+events_all <- left_join(events_all, select(Baseline, gupi, GOSE6mo, GOSE6mo.derived, GOS, ICPDevice, DecompressiveCran), by="gupi")
+
+
+Corrplots <- function(Corrtemp, title, alfa)
+{
+  corr.n_plot <- ggplot(Corrtemp, aes(x=icp.threshold, y=time.threshold, z=corr.n, fill=corr.n)) +
+    geom_raster(interpolate=TRUE, alpha = alfa) +
+    geom_contour(color="grey", alpha=0.5, size=0.5) +
+    stat_contour(breaks=c(0), size=1, colour="black") +
+    scale_fill_gradientn(colours = rev(jet.colors(100)),limits=c(-1,1),name= "Correlation", guide = guide_colorbar(reverse=TRUE)) +
+    theme_minimal() +
+    theme(legend.title = element_text(face="italic"))+#, text = element_text(size = 14)) +
+    xlab("Intensity greater than (mmHg)") +
+    ylab("Duration greater than (mins)") +
+    ggtitle(title)
+  
+  
+  return(corr.n_plot)
+}
+
+# Generate plot with SD:
+Corrplots2 <- function(Corrtemp, title, alfa, sdwidth)
+{
+  Corrtemp <- Corrtemp %>% mutate(tl.high=corr.n+sdwidth*sd.corr,
+                                  tl.low=corr.n-sdwidth*sd.corr)
+  
+  #Corrtemp <- Corrtemp %>% mutate(tl.high = ifelse()
+  corr.n_plot <- ggplot(Corrtemp, aes(x=icp.threshold, y=time.threshold, z=corr.n, fill=corr.n)) +
+    geom_raster(interpolate=TRUE, alpha = alfa) +
+    geom_contour(color="grey", alpha=0.5, size=0.5) +
+    stat_contour(breaks=c(0), size=1, colour="black") +
+    stat_contour(breaks = c(0), size=1, aes(x=icp.threshold, y=time.threshold, z=tl.low), colour="grey") +
+    stat_contour(breaks = c(0), size=1, aes(x=icp.threshold, y=time.threshold, z=tl.high), colour="white") +
+    scale_fill_gradientn(colours = rev(jet.colors(100)),limits=c(-1,1),name= "Correlation", guide = guide_colorbar(reverse=TRUE)) +
+    theme_minimal() +
+    theme(legend.title = element_text(face="italic"))+#, text = element_text(size = 14)) +
+    xlab("Intensity greater than (mmHg)") +
+    ylab("Duration greater than (mins)") +
+    ggtitle(paste0(title, ", ", "+-", sdwidth, " SD"))
+  
+  
+  return(corr.n_plot)
+}
+
+# -----------------------------------------------------
+
+
+
+# colors <- scale_color_manual(values=brewer.pal(n=8, name="Spectral"))
+
+
+Thresholds <- c(0,10,15,20,25,30)
+
+
+
+
+# dev.new()
+# print(ggplot(PTD$tot$PTD.long %>% filter(Threshold %in% Thresholds), aes(x=mmHgH)) + geom_histogram() + facet_wrap(~Threshold, scales="free_x"))
+
+g <- "C:/Users/Cecilia KI/Dropbox/Cecilia/ICP burden/Manuscript/PLOS submission/Revision/Figures/"
+
+
+# Figure 1: ICP profile of one patient:
+
+temp <- ms %>% filter(gupi == unique(ms$gupi)[9]) %>% mutate(Days = difftime(datetime, "1970-01-01 00:00", units="days"))
+
+tiff(file=paste0(g,"Fig1.tiff"), width=5.2, height=3,units='in', res=300)
+
+print(ggplot(temp, aes(x=Days, y=icp, y2=10)) +
+        geom_flame(fill = "deepskyblue", alpha=0.5) +
+        geom_line(size=0.1) + geom_hline(yintercept = 10, size=0.5) +
+        ylim(0,max(temp$icp)) +
+        xlim(2,max(temp$Days)) +
+        theme_minimal() +
+        theme(legend.title = element_text(face="italic")) +
+        xlab("Days post injury") + ylab("ICP"))# +
+        #text = element_text(size = 14), axis.text = element_text(size=14))
+        #ggtitle("ICP monitoring from one patient") + 
+
+dev.off()
+
+
+# Fig 2: GOSE distributions:
+tiff(file=paste0(g,"Fig2.tiff"), width=2.63, height=2.63,units='in', res=300)
+
+print(ggplot(Baseline, aes(GOSE6mo)) + geom_bar(fill = "darkgrey") + xlab("GOS-E score at 6 months") + ylab("Count") +
+        scale_x_discrete(limits=c(1:8)) +
+        theme_minimal()) #+
+        #theme(legend.title = element_text(face="italic")), text = element_text(size = 14))
+dev.off()
+
+
+# Fig 3:
+tiff(file=paste0(g,"Fig3.tiff"), width=7.5, height=4,units='in', res=300)
+
+ggarrange(Corrplots(Corr.list$ICP$All, paste0("A. All patients, n=", nrow(Baseline)), 1),
+             Corrplots2(Corr.list$ICP$Bootstraps$All, "B. 1000 bootstraps", 1, 2), nrow=1, common.legend = TRUE, legend="right")
+dev.off()
+
+
+
+# Fig 4:
+tiff(file=paste0(g,"Fig4.tiff"), width=7.5, height=4,units='in', res=300)
+#dev.new()
+ggarrange(Corrplots(Corr.list$ICP$PRxBelow0.3, "A. Intact autoregulation", 1),
+             Corrplots(Corr.list$ICP$PRxAbove0.3, "B. Impaired autoregulation", 1), nrow=1, common.legend = TRUE, legend="right")
+dev.off()
+
+
+# Fig 5:
+tiff(file=paste0(g,"Fig5.tiff"), width=2.63, height=2.63,units='in', res=300)
+
+print(ggplot(Autoregtimes, aes(x=n, y=Activetime/Monitortime)) + geom_line() +
+        theme_minimal() +
+        theme(legend.title = element_text(face="italic")) +
+        xlim(0,200) +
+        xlab("Patient number") + ylab("Time with intact autoregulation\n /Total monitoring time"))
+dev.off()
+
+
+# Fig 6, PTD stratified per GOSE:
+tiff(file=paste0(g,"Fig6.tiff"), width=5.2, height=3,units='in', res=300)
+
+print(ggplot(data=PTD$tot$PTD.GOSE, aes(Threshold, mean.mmHgH, group=GOSE6mo, color=GOSE6mo)) + geom_line() +
+        scale_color_manual(values=rev(jet.colors(8)), name="6 month GOS-E") +
+        #scale_color_brewer(palette="RdYlBu") +
+        theme_minimal() +
+        theme(legend.title = element_text(face="italic")) +#, text=element_text(size=12), axis.text=element_text(size = 14)) +
+        #xlim(0,30) +
+        #ggtitle("Mean PTD stratified by 6-month outcome") +
+        xlab("ICP threshold (mmHg)") + ylab("PTD (mmHg·h), mean"))
+
+dev.off()
+
+
+# Fig 7: Mean PTD, favourable/unfavourable outcome:
+
+PTD.fav.plot <- ggplot(PTD$tot$PTD.FAV, aes(x=Threshold, y=mean.mmHgH, color=Outcome)) + geom_line() +
+  scale_color_manual(values=c("deepskyblue", "orange3"), name="Outcome") +
+  # scale_color_brewer(palette="Accent", name="Outcome") +
+  theme_minimal() +
+  theme(legend.title = element_text(face="italic")) +
+  ggtitle("A. Favourable vs Unfavourable outcome") + xlab("ICP threshold") + ylab("PTD (Mean mmHg·h)")
+
+PTD.mort.plot <- ggplot(PTD$tot$PTD.Mort, aes(x=Threshold, y=mean.mmHgH, color=Mort6mo)) + geom_line() +
+  scale_color_manual(values=c("darkgreen", "red"), name="6 month mortality") +
+  #labs(color = "6 month mortality") +
+  theme_minimal() +
+  theme(legend.title = element_text(face="italic")) +
+  ggtitle("B. 6 month mortality") + xlab("ICP threshold") + ylab("PTD (Mean mmHg·h)")
+
+tiff(file=paste0(g,"Fig7.tiff"), width=7.5, height=3,units='in', res=300)
+ggarrange(PTD.fav.plot, PTD.mort.plot, ncol=2)#, common.legend = TRUE, legend="right")
+dev.off()
+
+
+# Fig 8:
+Active.Fav.plot <- ggplot(PTD$active$PTD.FAV, aes(x=Threshold, y=mean.mmHgH, color=Outcome)) +
+  geom_line() +
+  scale_color_manual(values=c("deepskyblue", "orange3"), name="Outcome") +
+  theme_minimal() +
+  theme(legend.title = element_text(face="italic")) +
+  ggtitle("A. Intact autoregulation") + xlab("ICP threshold") + ylab("PTD (Mean mmHg·h)")
+
+Passive.Fav.plot <- ggplot(PTD$passive$PTD.FAV, aes(x=Threshold, y=mean.mmHgH, color=Outcome)) +
+  geom_line() +
+  scale_color_manual(values=c("deepskyblue", "orange3"), name="Outcome") +
+  theme_minimal() +
+  theme(legend.title = element_text(face="italic")) +
+  ggtitle("B. Impaired autoregulation") + xlab("ICP threshold") + ylab("PTD (Mean mmHg·h)")
+
+Active.Mort.plot <- ggplot(PTD$active$PTD.Mort, aes(x=Threshold, y=mean.mmHgH, color=Mort6mo)) +
+  geom_line() +
+  scale_color_manual(values=c("darkgreen", "red"), name="Mortality") +
+  theme_minimal() +
+  theme(legend.title = element_text(face="italic")) +
+  ggtitle("C. Intact autoregulation") + xlab("ICP threshold") + ylab("PTD (Mean mmHg·h)")
+
+Passive.Mort.plot <- ggplot(PTD$passive$PTD.Mort, aes(x=Threshold, y=mean.mmHgH, color=Mort6mo)) +
+  geom_line() +
+  scale_color_manual(values=c("darkgreen", "red"), name="Mortality") +
+  theme_minimal() +
+  theme(legend.title = element_text(face="italic")) +
+  ggtitle("D. Impaired autoregulation") + xlab("ICP threshold") + ylab("PTD (Mean mmHg·h)")
+
+
+tiff(file=paste0(g,"Fig8.tiff"), width=6, height=6, units='in',res=300)
+
+ggarrange(annotate_figure(ggarrange(Active.Fav.plot, Passive.Fav.plot, ncol=2, common.legend = TRUE, legend="right"), top=text_grob("Favourable vs Unfavourable outcome", face="bold")),
+             annotate_figure(ggarrange(Active.Mort.plot, Passive.Mort.plot, ncol=2, common.legend = TRUE, legend="right"), top=text_grob("6 month mortality", face="bold")), nrow=2)
+dev.off()
+
+
+
+
+
+# Supplement Fig 1:
+tiff(file=paste0(g,"S1_Fig.tiff"), width=7.5, height=3,units='in', res=300)
+print(Corrplots(Corr_bootstraps.GOSE, "Ten randomly selected bootstraps", 1) + facet_wrap(~bootstrap, nrow=2)) 
+dev.off()
+
+
+# Supplement Fig 2:
+NoEVDDC <- length(unique((events_all %>% filter(!gupi %in% unique(filter(ms, !is.na(icp.evd))$gupi), gupi %in% filter(Baseline, DecompressiveCran==0)$gupi))$gupi))
+
+tiff(file=paste0(g,"S2_Fig.tiff"), width=7.5, height=7.5, units='in',res=300)
+
+ggarrange(Corrplots(Corr.list$ICP$All, paste0("A. All patients, n=",nrow(Baseline)), 1),
+         Corrplots(Corr.list$ICP$NoDC$All, paste0("B. No DC, n=", nrow(Baseline %>% filter(DecompressiveCran==0))), 1),
+         Corrplots(Corr.list$ICP$NoEVD$All, paste0("C. No EVD, n=", length(unique(ms$gupi))-length(unique(filter(ms, !is.na(icp.evd))$gupi))), 1),
+         Corrplots(Corr.list$ICP$NoDCNoEVD$All, paste0("D. No EVD, no DC, n=",NoEVDDC ), 1), nrow=2, ncol=2, common.legend=TRUE, legend="right")
+
+dev.off()
+
+
+# Supplement Fig 3, GOS:
+tiff(file=paste0(g,"S3_Fig.tiff"), width=7.5, height=7.5, units='in',res=300)
+
+ggarrange(Corrplots(Corr.list$GOS$ICP$All %>% filter(icp.threshold>=10, time.threshold>=5), title=paste0("A. GOS, n=", length(unique(Baseline$gupi))), alfa=1),
+          Corrplots(Corr.list$GOS$ICP$PRxBelow0.3 %>% filter(icp.threshold>=10, time.threshold>=5), title="B. GOS, Active autoregulation", alfa=1),
+          Corrplots(Corr.list$GOS$ICP$PRxAbove0.3 %>% filter(icp.threshold>=10, time.threshold>=5), title="C. GOS, Passive autoregulation", alfa=1),
+          nrow=2, ncol=2, common.legend = TRUE, legend="right")
+dev.off()
+
+
+# Suppl Fig 4, distr of monitoring time:
+tiff(file=paste0(g,"S4_Fig.tiff"), width=5.2, height=4, units='in',res=300)
+
+distr <- ggplot(data=Baseline, aes(x=monitor.time, fill=Mort6mo)) +
+  geom_histogram(alpha=0.5) +
+  scale_fill_manual(values=c("darkgreen", "red"), name="6 Month Mortality") +
+  scale_x_continuous(name="Monitoring time (days)") +
+  theme_minimal() + 
+  theme(legend.title = element_text(face="italic"))
+
+print(distr)
+dev.off()
+
+# Supplement Fig 5, Distribution of PTD:
+tiff(file=paste0(g,"S5_Fig.tiff"), width=7.5, height=5, units='in',res=300)
+distr <- ggplot(data=PTD$tot$PTD.long %>% filter(Threshold %in% Thresholds), aes(x=mmHgH, fill=as.factor(Mort6mo))) +
+  geom_histogram(alpha=0.5) +
+  scale_fill_manual(values=c("darkgreen", "red"), name="6 month mortality") +
+  facet_wrap(~Threshold, scales = "free") +
+  theme_minimal() +
+  theme(legend.title = element_text(face="italic")) +
+  xlab("PTD (mmHg·h)")
+
+print(distr)
+dev.off()
+