Upload New File

ICP_burden_201006.R

+rm(list=ls())
+
+library(dplyr)
+library(tidyr)
+library(ggplot2)
+library(gridExtra)
+
+#For parallel processing:
+library(foreach)
+library(doParallel)
+
+
+nas <- c("", NA)
+f <- "C:/Users/Cecilia KI/Dropbox/Cecilia/ICP burden/ICP_tempdata/"
+
+# New data: ---------------------
+Baseline <-  read.csv("C:/Users/Cecilia KI/Dropbox/Cecilia/ICP burden/Baseline.new_200625.csv", stringsAsFactors = FALSE) #Created in script ICP_burden_preparation_200625.R
+
+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),
+                    GOSE6mo.derived = as.factor(GOSE6mo.derived))
+
+ms <- ms %>% mutate(icp= ifelse(is.na(icp), icp.evd, icp))
+
+jet.colors <- colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))
+cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
+
+pt <- 10:40
+tth <- 5:360
+
+
+# Functions: --------------------------------
+# ICP_events_all:
+
+
+ICP_events_all <- function(x, events, pt, tth)
+{
+  ndays <- as.numeric(difftime(max(x$datetime), min(x$datetime), units="days"))
+  
+  
+  for (i in pt)
+  {
+    print(paste0("Patient ", j, ", i ",i))
+    x <- mutate(x, icp_up=ifelse(icp<i, 0, 1), index=1:nrow(x), icp_trans=c(NA, diff(icp_up)))
+    x[1,"icp_trans"] <- ifelse(x[1,"icp"]<i,0,1)
+    
+    #Get start and end times of events:
+    x_start <- filter(x, icp_trans==1) %>% select(datetime) %>% rename(time_start=datetime) #get starttimes of events
+    
+    x3 <- filter(x, icp_trans==-1)
+    endtime <- x3$index-1
+    x_end <- filter(x, index %in% endtime) %>% select(datetime) %>% rename(time_stop=datetime) #get endtime of events
+    
+    #make x_start, x_end same length:
+    if(nrow(x_end)<nrow(x_start))
+    {x_end <- rbind(x_end, NA)}
+    
+    # # Lines added 191207:
+    # if(is.na(x_end[nrow(x_end)]))
+    # {x_end[nrow(x_end)] <- max(x$datetime)}
+    
+    # if(nrow(x_start)<nrow(x_end))
+    # {x_start <- rbind(NA, x_start)}
+    
+    
+    h <- cbind(x_start, x_end)
+    colnames(h) <- c("time_start", "time_stop") #To make consistent colnames if no start or end times
+    h <- h %>% mutate(PRx.mean = NA, PAx.mean = NA, RAC.mean = NA, RAP.mean = NA)
+    
+    #Get length of each event:
+    h$timediff <- difftime(h$time_stop, h$time_start, units="mins") #length of each ICP event
+    
+    if(nrow(h) != 0)
+    {
+      for (k in 1:nrow(h))
+        {
+        h$PRx.mean[k] <- mean((x %>% filter(datetime >= h$time_start[k], datetime <= h$time_stop[k]))$PRx, na.rm = TRUE)
+        h$PAx.mean[k] <- mean((x %>% filter(datetime >= h$time_start[k], datetime <= h$time_stop[k]))$PAx, na.rm = TRUE)
+        h$RAC.mean[k] <- mean((x %>% filter(datetime >= h$time_start[k], datetime <= h$time_stop[k]))$RAC, na.rm = TRUE)
+        h$RAP.mean[k] <- mean((x %>% filter(datetime >= h$time_start[k], datetime <= h$time_stop[k]))$RAP, na.rm = TRUE)
+      }
+    }
+    
+    events.temp <- data.frame(PatientID=j, icp.threshold=rep(i, length(tth)), time.threshold=tth, n=NA, n.fixed=NA, n.fixed.exact=NA)
+    
+    for (t in 1:nrow(events.temp))
+    {
+      print(paste0("t",t))
+      tt <- events.temp$time.threshold[t]
+      h.temp <- h %>% filter(timediff >= tt) %>%
+                      mutate(n.fixed = floor(timediff/tt),
+                             n.fixed.exact = timediff/tt)
+      
+      events.temp$n[t] <- nrow(h.temp)
+      events.temp$n.fixed[t] <- sum(floor(h.temp$n.fixed))
+      events.temp$n.fixed.exact[t] <- sum(h.temp$n.fixed.exact)
+      events.temp$n.PRx.active[t] <- nrow(filter(h.temp, PRx.mean <= 0))
+      events.temp$n.PRx.passive[t] <- nrow(filter(h.temp, PRx.mean > 0))
+      events.temp$n.PAx.active[t] <- nrow(filter(h.temp, PAx.mean <= 0))
+      events.temp$n.PAx.passive[t] <- nrow(filter(h.temp, PAx.mean > 0))
+      events.temp$n.RAC.active[t] <- nrow(filter(h.temp, RAC.mean <= 0))
+      events.temp$n.RAC.passive[t] <- nrow(filter(h.temp, RAC.mean > 0))
+      events.temp$n.RAP.active[t] <- nrow(filter(h.temp, RAP.mean <= 0))
+      events.temp$n.RAP.passive[t] <- nrow(filter(h.temp, RAP.mean > 0))
+      events.temp$PRx.above.0.3[t] <- nrow(filter(h.temp, PRx.mean > 0.3))
+      events.temp$PRx.below.0.3[t] <- nrow(filter(h.temp, PRx.mean <= 0.3))
+      events.temp$PRx.above.0.2[t] <- nrow(filter(h.temp, PRx.mean > 0.2))
+      events.temp$PRx.below.0.2[t] <- nrow(filter(h.temp, PRx.mean <= 0.2))
+      events.temp$PRx.above.0.1[t] <- nrow(filter(h.temp, PRx.mean > 0.1))
+      events.temp$PRx.below.0.1[t] <- nrow(filter(h.temp, PRx.mean <= 0.1))
+      events.temp$PRx.mean[t] <- sum(h.temp$PRx.mean*as.integer(h.temp$timediff))/sum(as.integer(h.temp$timediff))
+      events.temp$Mean.PRx.mean[t] <- mean(h.temp$PRx.mean)
+      }
+    
+    events.temp <- events.temp %>% mutate(days = ndays, 
+                                          n.per.day=n/days,
+                                          n.fixed.per.day=n.fixed/days,
+                                          n.fixed.exact.per.day=n.fixed.exact/days)
+    events <- rbind(events, events.temp)
+  }
+  return(events)
+}
+
+# events <- temp %>% filter(TIL.Tier==0) %>% rename(Outcome=GOSE6mo)
+# m="n"
+
+ICP_correlations <- function(events, m)
+{
+  events <- events %>% select(icp.threshold, time.threshold, Outcome, m)
+  events <- events %>% rename(k = m)
+  events.summary <- events %>% group_by(icp.threshold, time.threshold, Outcome) %>%
+                               summarise(mean.n=mean(k, na.rm=TRUE)) %>%
+                               filter(!is.na(Outcome)) # %>%
+                               #ungroup(events)
+  
+  # new lines: 191210:
+  temp <- data.frame(Outcome = rep(unique(events$Outcome), each=length(pt)), 
+                     icp.threshold = pt)
+  temp2 <- data.frame(icp.threshold=rep(pt, each=length(tth), ncol=length(tth)))
+  temp2 <- cbind(temp2, time.threshold = tth)
+  temp <- left_join(temp, temp2, by="icp.threshold")
+  
+  temp <- left_join(temp, events.summary, by=c("icp.threshold", "time.threshold", "Outcome"))
+  
+  temp[is.na(temp)] <- 0
+  events.summary <- temp
+  
+  # line added 200128:
+  events.summary$Outcome <- as.numeric(events.summary$Outcome)
+  
+  #Calculate correlation coefficients between outcome and number of events:
+  Corr.temp <- events.summary %>% group_by(icp.threshold, time.threshold) %>%
+               summarise(corr.n = as.numeric(cor.test(Outcome, mean.n)$estimate),
+                         min.CI = min(as.numeric(cor.test(Outcome, mean.n, method = "pearson")$conf.int)),
+                         max.CI = max(as.numeric(cor.test(Outcome, mean.n, method = "pearson")$conf.int)))
+ 
+  n.patients <- events %>% filter(k != 0) %>% group_by(icp.threshold, time.threshold) %>% summarise(n.pat.nonzero = n()) #%>% ungroup(events)
+  
+  Corr.temp <- left_join(Corr.temp, n.patients, by=c("icp.threshold", "time.threshold"))
+  
+
+  return(Corr.temp)
+}
+
+# Function to generate plots:
+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=2, 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=2, colour="black") +
+    stat_contour(breaks = c(0), size=2, aes(x=icp.threshold, y=time.threshold, z=tl.low), colour="grey") +
+    stat_contour(breaks = c(0), size=2, 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)
+}
+
+# Calculate events and Bootstrappings and correlations for bootstrapping: ---------------------
+
+# Create events (ms_200625 was used): ---------------------
+
+# # Create ICP events: -------------------------------
+#Thresholds:
+events <- data.frame()
+Patients <- unique(Baseline$gupi)
+
+# Create events list, filtered data:
+no_cores <- detectCores()
+clust <- makeCluster(no_cores-1)
+registerDoParallel(clust)
+
+
+StartTime <- Sys.time()
+events_all <- foreach(j = 1:length(Patients), .combine = rbind, .packages = "dplyr") %dopar% {ICP_events_all(filter(ms, gupi == Patients[j]), events, pt, tth)}
+Sys.time() - StartTime
+stopCluster(clust)
+
+events_all <- events_all %>% mutate(gupi = Patients[PatientID]) %>% select(gupi, everything(), -PatientID)
+# write.csv(events_all, "C:/Users/Cecilia KI/Dropbox/Cecilia/ICP burden/ICP_tempdata/200625/events_all_200625.csv", row.names = FALSE)
+
+
+# Bootstrap correlations, 1000 times: ------------------------------------------------------
+
+# Bootstrap events, all patients:
+load(paste0(f,"/200625/bootstrapping.RData")) # Created in ICP_burden_bootstraps_200625.R (On server)
+
+# -----------------------------------
+
+
+# Make list of correlations: ----------------------------------------------
+events_all <- read.csv(paste0(f,"200625/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")
+
+Corr.list <- list()
+Corr.list$ICP$All <- ICP_correlations(events_all %>% rename(Outcome = GOSE6mo), m="n")
+Corr.list$ICP$Bootstraps$All <- read.csv(paste0(f,"200625/Bootstraps.GOSE6mo.summary.csv"), stringsAsFactors = FALSE, na.strings = nas)
+Corr.list$ICP$Bootstraps$NoDC <- read.csv(paste0(f,"200625/Bootstraps.NoDC.GOSE6mo.summary.csv"), stringsAsFactors = FALSE, na.strings = nas)
+Corr.list$ICP$Bootstraps$DC <- read.csv(paste0(f,"200625/Bootstraps.DC.GOSE6mo.summary.csv"), stringsAsFactors = FALSE, na.strings = nas)
+Corr.list$ICP$Bootstraps$EVD <- read.csv(paste0(f,"200625/Bootstraps.EVD.GOSE6mo.summary.csv"), stringsAsFactors = FALSE, na.strings = nas)
+
+Corr.list$ICP$PRxBelow0.3 <- ICP_correlations(events_all %>% rename(Outcome = GOSE6mo), "PRx.below.0.3")
+Corr.list$ICP$PRxAbove0.3 <- ICP_correlations(events_all %>% rename(Outcome = GOSE6mo), "PRx.above.0.3")
+Corr.list$ICP$DC$All <- ICP_correlations(events_all %>% filter(DecompressiveCran==1) %>% rename(Outcome=GOSE6mo), m="n")
+Corr.list$ICP$NoDC$All <- ICP_correlations(events_all %>% filter(DecompressiveCran==0) %>% rename(Outcome=GOSE6mo), m="n")
+Corr.list$ICP$NoDC$PRxBelow0.3 <- ICP_correlations(events_all %>% filter(DecompressiveCran==0) %>% rename(Outcome = GOSE6mo), "PRx.below.0.3")
+Corr.list$ICP$NoDC$PRxAbove0.3 <- ICP_correlations(events_all %>% filter(DecompressiveCran==0) %>% rename(Outcome = GOSE6mo), "PRx.above.0.3")
+Corr.list$ICP$EVD$All <- ICP_correlations(events_all %>% filter(gupi %in% unique(filter(ms, !is.na(icp.evd))$gupi)) %>% rename(Outcome = GOSE6mo, m="n"))
+Corr.list$ICP$NoEVD$All <- ICP_correlations(events_all %>% filter(!gupi %in% unique(filter(ms, !is.na(icp.evd))$gupi)) %>% rename(Outcome = GOSE6mo, m="n"))
+
+Corr.list$ICP$NoDCNoEVD$All <- ICP_correlations(events_all %>% filter(!gupi %in% unique(filter(ms, !is.na(icp.evd))$gupi), gupi %in% filter(Baseline, DecompressiveCran==0)$gupi) %>% rename(Outcome=GOSE6mo), m="n")
+
+Corr.list$GOS$ICP$All <- ICP_correlations(events_all %>% rename(Outcome = GOS), m="n")
+Corr.list$GOS$ICP$PRxBelow0.3 <- ICP_correlations(events_all %>% rename(Outcome = GOS), m="PRx.below.0.3")
+Corr.list$GOS$ICP$PRxAbove0.3 <- ICP_correlations(events_all %>% rename(Outcome = GOS), m="PRx.above.0.3")
+
+View(Corr.list)
+# save(Corr.list, file = paste0(f,"200625/Corr.list.RData"))
+
+# --------------------------------------------------------------------------------------------------------------------------------------------------------
+load(paste0(f,"200625/Corr.list.RData"))
+Corr_bootstraps.GOSE <- read.csv(paste0(f,"200625/Corr_bootstraps.GOSE6mo.df.10.csv"), stringsAsFactors = FALSE)
+
+# Figures 200625:
+
+# Figure 1: ICP profile of one patient:
+library(heatwaveR)
+
+temp <- ms %>% filter(gupi == unique(ms$gupi)[9]) %>% mutate(Days = difftime(datetime, "1970-01-01 00:00", units="days"))
+dev.new()
+print(ggplot(temp, aes(x=Days, y=icp, y2=15)) +
+        geom_flame(fill = "grey") +
+        geom_line() + geom_hline(yintercept = 15, size=1) +
+        ylim(0,max(temp$icp)) +
+        xlim(2,max(temp$Days)) +
+        theme_minimal() +
+        theme(legend.title = element_text(face="italic"), text = element_text(size = 14), axis.text = element_text(size=14)) +
+        ggtitle("ICP monitoring from one patient") + xlab("Days post injury") + ylab("ICP"))
+
+
+
+# Fig 2: GOSE distributions:
+dev.new()
+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)))
+
+
+# Fig 3:
+dev.new()
+grid.arrange(Corrplots(Corr.list$ICP$All, paste0("A. All patients, n=", nrow(Baseline)), 1),
+             Corrplots2(Corr.list$ICP$Bootstraps$All, "B. 1000 bootstraps, n=227", 1, 2), nrow=1)
+
+# Fig 4:
+
+dev.new()
+print(Corrplots(Corr_bootstraps.GOSE, "Ten randomly selected bootstraps", 1) + facet_wrap(~bootstrap, nrow=2)) 
+
+# Fig 6:
+dev.new()
+grid.arrange(Corrplots(Corr.list$ICP$PRxBelow0.3, "A. Intact autoregulation", 1),
+             Corrplots(Corr.list$ICP$PRxAbove0.3, "B. Impaired autoregulation", 1), nrow=1)
+
+
+# Fig 5, active/passive time:
+# Autoregulation % active/passive:
+# % impaired autoregulation
+nrow(ms %>% filter(PRx.evd > 0.3 | PRx > 0.3))/nrow(ms)
+
+Autoregtimes <- data.frame(gupi=character(), Activetime = numeric(), Monitortime = numeric())
+
+for (i in 1:length(unique(Baseline$gupi)))
+{
+  print(Baseline$gupi[i])
+  x <- ms %>% filter(gupi == Baseline$gupi[i], !is.na(PRx)) %>% select(gupi, datetime, PRx)
+  
+  if(nrow(x) == 0)
+  {
+    activetime <- 0
+    monitortime <- 0
+  }
+  
+  else{
+  
+  x <- mutate(x, PRx.active=ifelse(PRx<0.3, 1, 0), index=1:nrow(x), prx_trans=c(NA, diff(PRx.active)))
+  x[1,"prx_trans"] <- ifelse(x[1,"PRx"]<0.3,1,-1)
+  
+  #Get start and end times of events, active:
+  x_start <- filter(x, prx_trans==1) %>% select(datetime) %>% rename(time_start=datetime) #get starttimes of events
+  x_start <- ungroup(x_start) %>% select(-gupi)
+  
+  x3 <- filter(x, prx_trans==-1)
+  endtime <- x3$index-1
+  x_end <- filter(x, index %in% endtime) %>% select(datetime) %>% rename(time_stop=datetime) #get endtime of events
+  x_end <- ungroup(x_end) %>% select(-gupi)
+  
+  #make x_start, x_end same length:
+  if(nrow(x_end)<nrow(x_start))
+  {x_end <- rbind(x_end, NA)}
+  
+  h <- cbind(x_start, x_end)
+  colnames(h) <- c("time_start", "time_stop") #To make consistent colnames if no start or end times
+  
+  if(is.na(h$time_stop[nrow(h)]))
+  {h$time_stop[nrow(h)] <- max(x$datetime)}
+  #Get time with active prx:
+  h$activetime <- difftime(h$time_stop, h$time_start, units="mins") #length of each ICP event
+  
+  activetime <- as.numeric(sum(h$activetime))
+  monitortime <- as.numeric(difftime(max(x$datetime), min(x$datetime), units = "mins"))
+  }
+  
+  Autoregtimes <- rbind(Autoregtimes, data.frame(gupi = Baseline$gupi[i], Activetime = activetime, Monitortime = monitortime))
+  }
+
+
+Autoregtimes <- arrange(Autoregtimes, Activetime/Monitortime)
+Autoregtimes$gupi <- factor(Autoregtimes$gupi, levels = Autoregtimes$gupi)
+Autoregtimes$n <- c(1:nrow(Autoregtimes))
+# write.csv(Autoregtimes, paste0(f,"200625/Autoregtimes.csv"), row.names=FALSE)
+
+# dev.new()
+# print(ggplot(Autoregtimes) + geom_bar(aes(x=as.factor(gupi), y=Monitortime/Monitortime), fill="darkgreen", stat="identity") +
+#         geom_bar(aes(x=as.factor(gupi), y=Activetime/Monitortime), fill="green", stat="identity")
+#       + ggtitle("Dark green = passive autoregulation, light green = active"))
+
+dev.new()
+print(ggplot(Autoregtimes, aes(x=n, y=Activetime/Monitortime)) + geom_line() +
+        theme_minimal() +
+        theme(legend.title = element_text(face="italic"), text = element_text(size = 14), axis.text = element_text(size=14)) +
+        xlim(0,200) +
+        ggtitle("Intact cerebral autoregulation") + xlab("Patient number") + ylab("Time with intact autoregulation / Total monitoring time"))
+
+
+# To supplement:
+NoEVDDC <- length(unique((events_all %>% filter(!gupi %in% unique(filter(ms, !is.na(icp.evd))$gupi), gupi %in% filter(Baseline, DecompressiveCran==0)$gupi))$gupi))
+
+dev.new()
+grid.arrange(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. Other monitor than 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)
+
+
+# To supplement, GOS:
+
+dev.new()
+grid.arrange(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)
+
+# For reference: -----------------------------
+dev.new() 
+grid.arrange(Corrplots2(Corr.list$ICP$Bootstraps$EVD, paste0("C. EVD, 1000 bootstraps, n=", nrow(filter(Baseline, !is.na(icp.evd.mean)))),1, 2),
+             Corrplots2(Corr.list$ICP$Bootstraps$DC, paste0("C. DC, 1000 bootstraps, n=", nrow(filter(Baseline, DecompressiveCran==1))),1, 2), nrow=1)
+
+# TIL:
+temp <- left_join(events_all, select(Baseline, gupi, TIL.Tier, TIL.max), by="gupi")
+Corr.list2 <- list()
+Corr.list2$TIL$TILHigh <- ICP_correlations(temp %>% filter(TIL.max >= 11) %>% rename(Outcome=GOSE6mo), m="n")
+Corr.list2$TIL$TILLow <- ICP_correlations(temp %>% filter(TIL.max < 11) %>% rename(Outcome=GOSE6mo), m="n")
+
+dev.new()
+grid.arrange(Corrplots(Corr.list2$TIL$TILHigh, title=paste0("Max Daily TIL >= 11, n=", length(unique((temp %>% filter(TIL.max >= 11))$gupi))), alfa=1),
+             Corrplots(Corr.list2$TIL$TILLow, title=paste0("Max Daily TIL < 11, n=", length(unique((temp %>% filter(TIL.max < 11))$gupi))), alfa=1), nrow=1)
+
+# --------------------------------------------
+
+# Distribution of monitoring time: ----------
+Baseline <- Baseline %>% mutate(Mort6mo = ifelse(DeathDate <= "1970-06-31", 1, 0))
+Baseline <- Baseline %>% mutate(Mort6mo = as.factor(ifelse(is.na(Mort6mo), 0, Mort6mo)))
+
+distr <- ggplot(data=Baseline, aes(x=monitor.time, fill=Mort6mo)) +
+  geom_histogram(alpha=0.5) +
+  scale_fill_manual(values=c("deepskyblue", "orange3"), name="6 Month Mortality") +
+  scale_x_continuous(name="Monitoring time (days)") +
+  theme_minimal() + 
+  theme(legend.title = element_text(face="italic"), text = element_text(size = 14), axis.text = element_text(size=14)) +
+  ggtitle("Distribution of monitoring time")
+  
+
+dev.new()
+print(distr)
+
+# Determine cut-offs (transition lines) and time in red zone: --------------------------------------------------------------------------------------------
+
+alfa <- 1
+
+#View(ggplot_build(corr.n_plot)$data)
+
+Blackline <- Corrplots(Corrtemp =  Corr.list$ICP$Bootstraps$All, paste0("n= ", length(unique(Baseline$gupi))), 1)
+Greywhitelines <- Corrplots2(Corr.list$ICP$Bootstraps$All %>% select(-corr.n) %>% rename(corr.n = mean.corr),
+                             "Mean correlations, 1000 bootstraps", 1, sdwidth=2)
+
+View(ggplot_build(Blackline)$data)
+View(ggplot_build(Greywhitelines)$data[[5]])
+View(ggplot_build(Greywhitelines)$data[[3]])
+
+View(ggplot_build(Corrplots(Corrtemp = Corr.list$ICP$PRxBelow0.3, "x",1)))
+
+# Thresholds for transition line:
+# +2*sd:
+ub.tl <- ggplot_build(Greywhitelines)$data[[5]] %>% mutate(ints=cut(x, breaks=seq(4.5,40.5, 1))) %>% group_by(ints) %>% summarise(time.threshold=median(y))
+ub.tl <- ub.tl %>% mutate(icp.threshold=c(13:32)) %>% select(-ints)
+
+# For black transition line:
+black.tl <- ggplot_build(Greywhitelines)$data[[3]] %>% mutate(ints=cut(x,breaks=seq(4.5,40.5,1))) %>% group_by(ints) %>% summarise(time.threshold=median(y))
+black.tl <- black.tl %>% mutate(icp.threshold=c(10:18)) %>% select(-ints)
+
+# Time in red zone:
+ms2 <- data.frame()
+ms <- ms %>% mutate(red.zone = 0)
+
+ICP_events_all2 <- function(x, tl, ms2)
+{
+  print(x$gupi[1])
+  ndays <- as.numeric(difftime(max(x$datetime), min(x$datetime), units="days"))
+  
+  for (i in 1:nrow(tl))
+  {
+    #print(paste0("Patient ", j, ", i ",tl$icp.threshold[i]))
+    x <- mutate(x, icp_up=ifelse(icp<tl$icp.threshold[i], 0, 1), index=1:nrow(x), icp_trans=c(NA, diff(icp_up)))
+    x[1,"icp_trans"] <- ifelse(x[1,"icp"]<tl$icp.threshold[i],0,1)
+    
+    #Get start and end times of events:
+    x_start <- x %>% ungroup %>% filter(icp_trans==1) %>% select(datetime) %>% rename(time_start=datetime) #get starttimes of events
+    
+    x3 <- x %>% ungroup %>% filter(icp_trans==-1)
+    endtime <- x3$index-1
+    x_end <- x %>% ungroup %>% filter(index %in% endtime) %>% select(datetime) %>% rename(time_stop=datetime) #get endtime of events
+    
+    #make x_start, x_end same length:
+    if(nrow(x_end)<nrow(x_start))
+    {x_end <- rbind(x_end, NA)}
+    
+    h <- cbind(x_start, x_end)
+    
+    if (nrow(h) !=0)
+    {
+      colnames(h) <- c("time_start", "time_stop") #To make consistent colnames if no start or end times
+      # h <- h %>% mutate(PRx.mean = NA, PAx.mean = NA, RAC.mean = NA, RAP.mean = NA)
+      
+      if(is.na(h$time_stop[nrow(h)]))
+      {
+        h$time_stop[nrow(h)] <- max(x$datetime)
+        h$time_stop <- as.POSIXct(h$time_stop, origin="1970-01-01")
+      }
+      
+      #Get length of each event:
+      h$timediff <- difftime(h$time_stop, h$time_start, units="mins") #length of each ICP event
+      h <- h %>% mutate(time_start2 = time_start + tl$time.threshold[i]*60, timediff2 = timediff-tl$time.threshold[i])
+      h <- h %>% filter(timediff2 > 0)
+      
+      if(nrow(h) > 0)
+      {
+        for (k in 1:nrow(h))
+        {
+          x <- x %>% mutate(red.zone = ifelse(between(datetime, h$time_start2[k], h$time_stop[k]), 1,red.zone))
+        }
+      }
+      
+      
+     # print(nrow(x %>% filter(red.zone==1)))
+    }
+    
+    
+    
+    
+  }
+  
+  x <- x %>% select(-icp_up, -icp_trans, -index)
+  
+  ms2 <- rbind(ms2, x)
+  return(ms2)
+}
+
+for (j in 1:length(unique(ms$gupi)))
+{
+  x <- ms %>% filter(gupi==unique(ms$gupi)[j]) %>% select(gupi, datetime, icp, red.zone, GOSE6mo.derived)
+  ms2 <- ICP_events_all2(x, black.tl, ms2)
+}
+
+rz <- ms2 %>% group_by(gupi) %>% summarise(n.rz = sum(red.zone), rz=mean(red.zone))
+View(rz)
+Baseline <- left_join(Baseline, rz, by="gupi")
+View(Baseline)
+
+# write.csv(Baseline, "C:/Users/Cecilia KI/Dropbox/Cecilia/ICP burden/Baseline.new_200625.csv", row.names = FALSE) # red zone time added in ICP_burden_200625.R)
+
+
+# Time above upper 2 SD bound: ----------------------------------------
+ms3 <- data.frame()
+for(j in 1:length(unique(ms$gupi)))
+{
+  x <- ms %>% filter(gupi == unique(ms$gupi)[j]) %>% select(gupi, datetime, icp, red.zone, GOSE6mo.derived)
+  ms3 <- ICP_events_all2(x, ub.tl, ms3)
+}
+
+rz.ub <- ms3 %>% group_by(gupi) %>% summarise(n.rz.ub=sum(red.zone), rz.ub=mean(red.zone))
+
+Baseline <- left_join(Baseline, rz.ub, by="gupi")
+# write.csv(Baseline, "C:/Users/Cecilia KI/Dropbox/Cecilia/ICP burden/Baseline.new_200625.csv", row.names = FALSE) # red zone time added in ICP_burden_200625.R)
+
+# ----------------------------------------------
+
+