Upload New File

ICP_burden_AUC_201006.R

+
+# Script that calculates PTD for all patients:
+# Doses compared with Kolmogorov-Smirnov.
+# For PTD CPP, MAP, PRx, see ICP_burden_AUC_200127.R
+# For means, see ICP_burden_AUC_191130.csv (or 191205)
+
+
+rm(list=ls())
+
+library(ggplot2)
+library(gridExtra)
+library(dplyr)
+library(tidyr)
+
+library(rms)
+library(RColorBrewer)
+library(ROCR)
+
+nas = c(NA, "")
+f <- "C:/Users/Cecilia KI/Dropbox/Cecilia/ICP burden/ICP_tempdata/"
+
+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
+
+
+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 <- 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"))
+activetime <- ms %>% filter(PRx <= 0.3) %>% group_by(gupi) %>% summarise(at=n()/60)
+passivetime <- ms %>% filter(PRx > 0.3) %>% group_by(gupi) %>% summarise(pt=n()/60)
+Baseline <- left_join(Baseline, activetime, by="gupi")
+Baseline <- left_join(Baseline, passivetime, by="gupi")
+rm(activetime, passivetime)
+
+ms <- left_join(ms, Baseline %>% select(gupi, GOSE6mo.derived), by="gupi")
+ms <- ms %>% select(gupi, datetime, icp, CPP, PRx, GOSE6mo, GOS, TIL.mean, TIL.max, DecompressiveCran, DecompressiveCran.HRtime, monitor.time)
+
+
+# Calculate AUC for different thresholds of ICP: -------------------------
+AUC.list <- list()
+AUCfcn <- function(i, dta){dta %>% filter(icp >= i) %>% group_by(gupi) %>% summarise(AUC = round(sum(icp-i)/60, 2), 
+                                                                               AUC.hours = round(sum(icp-i)/60/monitor.time[1],2),
+                                                                               Hours =  n()/60)} # Bad naming
+
+AUC.list$tot <- lapply(0:40, AUCfcn, ms)
+AUC.list$active <- lapply(0:40, AUCfcn, dta = ms %>% filter(PRx < 0.3))
+AUC.list$passive <- lapply(0:40, AUCfcn, dta = ms %>% filter(PRx >= 0.3))
+
+AUC.list$active.normalized  <- lapply(0:40, AUCfcn, dta = ms %>% filter(PRx < 0.3))
+AUC.list$passive.normalized  <- lapply(0:40, AUCfcn, dta = ms %>% filter(PRx >= 0.3))
+
+
+for (i in 0:40)
+{
+  AUC.list$active.normalized[[i+1]] <- left_join(AUC.list$active.normalized[[i+1]], select(Baseline, gupi, at), by="gupi") %>% select(-AUC.hours) %>%
+                                       mutate(AUC.active.normhours = AUC/at) %>% select(-at)
+  
+  AUC.list$passive.normalized[[i+1]] <- left_join(AUC.list$passive.normalized[[i+1]], select(Baseline, gupi, pt), by="gupi") %>% select(-AUC.hours) %>%
+    mutate(AUC.active.normhours = AUC/pt) %>% select(-pt)
+}
+  
+# View(AUC.list$passive.normalized[[28]])  
+  
+for (i in 0:40)
+  {
+    colnames(AUC.list$tot[[i+1]]) <- c("gupi", paste0("PTD.ICP",i), paste0("Hours.",i), paste0("PTD.ICP",i,".perday"))
+    colnames(AUC.list$active[[i+1]]) <- c("gupi", paste0("PTD.ICP",i), paste0("Hours.",i), paste0("PTD.ICP",i,".perday"))
+    colnames(AUC.list$passive[[i+1]]) <- c("gupi", paste0("PTD.ICP",i), paste0("Hours.",i), paste0("PTD.ICP",i,".perday"))
+    colnames(AUC.list$active.normalized[[i+1]]) <- c("gupi", paste0("PTD.ICP",i), paste0("Hours.",i), paste0("PTD.ICP",i,".perhouractive"))
+    colnames(AUC.list$passive.normalized[[i+1]]) <- c("gupi", paste0("PTD.ICP",i), paste0("Hours.",i), paste0("PTD.ICP",i,".perhourpassive"))
+    
+  }
+
+
+
+# Make one data frame of all AUC:s:
+AUC.tot <- left_join(AUC.list$tot[[1]], AUC.list$tot[[2]], by="gupi")
+AUC.active <- left_join(AUC.list$active[[1]], AUC.list$active[[2]], by="gupi")
+AUC.passive <- left_join(AUC.list$passive[[1]], AUC.list$passive[[2]], by="gupi")
+AUC.active.norm <- left_join(AUC.list$active.normalized[[1]], AUC.list$active.normalized[[2]], by="gupi")
+AUC.passive.norm <- left_join(AUC.list$passive.normalized[[1]], AUC.list$passive.normalized[[2]], by="gupi")
+
+for (i in 3:length(AUC.list$tot))
+{
+  AUC.tot <- left_join(AUC.tot, AUC.list$tot[[i]], by="gupi")
+  AUC.active <- left_join(AUC.active, AUC.list$active[[i]], by="gupi")
+  AUC.passive <- left_join(AUC.passive, AUC.list$passive[[i]], by="gupi")
+  AUC.active.norm <- left_join(AUC.active.norm, AUC.list$active.normalized[[i]], by="gupi")
+  AUC.passive.norm <- left_join(AUC.passive.norm, AUC.list$passive.normalized[[i]], by="gupi")
+}
+
+
+AUC.tot[is.na(AUC.tot)] <- 0
+AUC.active[is.na(AUC.active)] <- 0
+AUC.passive[is.na(AUC.passive)] <- 0
+AUC.active.norm[is.na(AUC.active.norm)] <- 0
+AUC.passive.norm[is.na(AUC.passive.norm)] <- 0
+
+AUC.tot <- left_join(AUC.tot, select(Baseline, gupi, GOSE6mo, GOS, GOSE6mo.derived, Mort90d, Mort6mo, Outcome), by="gupi")
+AUC.active <- left_join(AUC.active, select(Baseline, gupi, GOSE6mo, GOS, GOSE6mo.derived, Mort90d, Mort6mo, Outcome), by="gupi")
+AUC.passive <- left_join(AUC.passive, select(Baseline, gupi, GOSE6mo, GOS, GOSE6mo.derived, Mort90d, Mort6mo, Outcome), by="gupi")
+AUC.active.norm <- left_join(AUC.active.norm, select(Baseline, gupi, at, GOSE6mo, GOS, GOSE6mo.derived, Mort90d, Mort6mo, Outcome), by="gupi")
+AUC.passive.norm <- left_join(AUC.passive.norm, select(Baseline, gupi, pt, GOSE6mo, GOS, GOSE6mo.derived, Mort90d, Mort6mo, Outcome), by="gupi")
+
+# View(left_join(AUC.passive.norm %>% select(gupi, PTD.ICP20:PTD.ICP23.perhourpassive), select(Baseline, gupi, at, pt), by="gupi") %>% mutate(at=round(at/60,2), pt=round(pt/60,2)))
+
+PTD.fcn <- function(AUC)
+{
+  PTDperday <- AUC %>% select(gupi, contains("day"), GOSE6mo)
+  PTD.perday.long <- gather(PTDperday, Threshold, mmHgH.perday, PTD.ICP0.perday:PTD.ICP40.perday) %>% mutate(Threshold = gsub("PTD.ICP","",Threshold))
+  PTD.perday.long$Threshold <- as.numeric(gsub(".perday", "", PTD.perday.long$Threshold))
+  
+  PTD <- AUC %>% select(-contains("day"))
+  PTD.long <- gather(PTD, Threshold, mmHgH, PTD.ICP0:PTD.ICP40) %>% mutate(Threshold = as.numeric(gsub("PTD.ICP","",Threshold)))
+  PTD.long <- left_join(PTD.perday.long, PTD.long)
+  
+  PTD.long <- PTD.long %>% mutate(mmHgH.perday = ifelse(is.na(mmHgH.perday), 0, mmHgH.perday),
+                                  mmHgH = ifelse(is.na(mmHgH), 0, mmHgH))
+  
+  PTD.summaries <- PTD.long %>% group_by(Threshold) %>% summarise(mean.mmHgH = mean(mmHgH, na.rm=TRUE), median.mmHgH = median(mmHgH, na.rm = TRUE),
+                                                                  mean.mmHgH.perday = mean(mmHgH.perday, na.rm=TRUE), median.mmHgH.perday = median(mmHgH.perday, na.rm=TRUE)) 
+  
+  PTD.GOSE <- PTD.long %>% group_by(GOSE6mo, Threshold) %>% summarise(mean.mmHgH = mean(mmHgH, na.rm=TRUE), median.mmHgH = median(mmHgH, na.rm = TRUE),
+                                                                   mean.mmHgH.perday = mean(mmHgH.perday, na.rm=TRUE), median.mmHgH.perday = median(mmHgH.perday, na.rm=TRUE)) 
+  
+  PTD.FAV <- PTD.long %>% group_by(Outcome, Threshold) %>% summarise(mean.mmHgH = mean(mmHgH, na.rm=TRUE), median.mmHgH = median(mmHgH, na.rm = TRUE),
+                                                                     mean.mmHgH.perday = mean(mmHgH.perday, na.rm=TRUE), median.mmHgH.perday=median(mmHgH.perday, na.rm=TRUE)) 
+  
+  PTD.Mort <- PTD.long %>% group_by(Mort6mo, Threshold) %>% summarise(mean.mmHgH = mean(mmHgH, na.rm=TRUE), median.mmHgH = median(mmHgH, na.rm = TRUE),
+                                                                      mean.mmHgH.perday = mean(mmHgH.perday, na.rm=TRUE), median.mmHgH.perday = median(mmHgH.perday, na.rm=TRUE)) 
+PTD.list <- list(PTD.long, PTD.summaries, PTD.GOSE, PTD.FAV, PTD.Mort)   
+names(PTD.list) <- c("PTD.long","PTD.summaries", "PTD.GOSE", "PTD.FAV", "PTD.Mort")
+return(PTD.list)  
+}
+
+PTD.fcn2 <- function(AUC)
+{
+  PTDpernormhour <- AUC %>% select(gupi, contains("perhour"), GOSE6mo)
+  colnames(PTDpernormhour) <- gsub("passive", "", colnames(PTDpernormhour))
+  colnames(PTDpernormhour) <- gsub("active", "", colnames(PTDpernormhour))
+  
+  
+  PTD.pernormhour.long <- gather(PTDpernormhour, Threshold, mmHgH.pernormhour, PTD.ICP0.perhour:PTD.ICP40.perhour) %>% mutate(Threshold = gsub("PTD.ICP","",Threshold))
+  PTD.pernormhour.long$Threshold <- as.numeric(gsub(".perhour", "", PTD.pernormhour.long$Threshold))
+  
+  PTD <- AUC %>% select(-contains("hour"))
+  PTD.long <- gather(PTD, Threshold, mmHgH, PTD.ICP0:PTD.ICP40) %>% mutate(Threshold = as.numeric(gsub("PTD.ICP","",Threshold)))
+  PTD.long <- left_join(PTD.pernormhour.long, PTD.long)
+  
+  PTD.long <- PTD.long %>% mutate(mmHgH.pernormhour = ifelse(is.na(mmHgH.pernormhour), 0, mmHgH.pernormhour),
+                                  mmHgH = ifelse(is.na(mmHgH), 0, mmHgH))
+  
+  PTD.summaries <- PTD.long %>% group_by(Threshold) %>% summarise(mean.mmHgH = mean(mmHgH, na.rm=TRUE), median.mmHgH = median(mmHgH, na.rm = TRUE),
+                                                                  mean.mmHgH.pernormhour = mean(mmHgH.pernormhour, na.rm=TRUE), median.mmHgH.pernormhour = median(mmHgH.pernormhour, na.rm=TRUE)) 
+  
+  PTD.GOSE <- PTD.long %>% group_by(GOSE6mo, Threshold) %>% summarise(mean.mmHgH = mean(mmHgH, na.rm=TRUE), median.mmHgH = median(mmHgH, na.rm = TRUE),
+                                                                      mean.mmHgH.pernormhour = mean(mmHgH.pernormhour, na.rm=TRUE), median.mmHgH.pernormhour = median(mmHgH.pernormhour, na.rm=TRUE)) 
+  
+  PTD.FAV <- PTD.long %>% group_by(Outcome, Threshold) %>% summarise(mean.mmHgH = mean(mmHgH, na.rm=TRUE), median.mmHgH = median(mmHgH, na.rm = TRUE),
+                                                                     mean.mmHgH.pernormhour = mean(mmHgH.pernormhour, na.rm=TRUE), median.mmHgH.pernormhour=median(mmHgH.pernormhour, na.rm=TRUE)) 
+  
+  PTD.Mort <- PTD.long %>% group_by(Mort6mo, Threshold) %>% summarise(mean.mmHgH = mean(mmHgH, na.rm=TRUE), median.mmHgH = median(mmHgH, na.rm = TRUE),
+                                                                      mean.mmHgH.pernormhour = mean(mmHgH.pernormhour, na.rm=TRUE), median.mmHgH.pernormhour = median(mmHgH.pernormhour, na.rm=TRUE)) 
+  PTD.list <- list(PTD.long, PTD.summaries, PTD.GOSE, PTD.FAV, PTD.Mort)   
+  names(PTD.list) <- c("PTD.long","PTD.summaries", "PTD.GOSE", "PTD.FAV", "PTD.Mort")
+  return(PTD.list)  
+}
+
+
+PTD <- list()
+PTD$tot <- PTD.fcn(AUC.tot)
+PTD$active <- PTD.fcn(AUC.active)
+PTD$passive <- PTD.fcn(AUC.passive)
+PTD$active.norm <- PTD.fcn2(AUC.active.norm)
+PTD$passive.norm <- PTD.fcn2(AUC.passive.norm)
+
+
+# save(PTD, file=paste0(f,"200625/PTD.RData"))
+# ----------------------------------------------------------------
+
+
+
+
+# Plots:-----------------------------------------------------------
+load(file=paste0(f,"200625/PTD.RData"))
+
+# colors <- scale_color_manual(values=brewer.pal(n=8, name="Spectral"))
+View(PTD)
+
+
+Thresholds <- c(0,10,15,20,25,30)
+
+# Distribution of PTD:
+distr <- ggplot(data=PTD$tot$PTD.long %>% filter(Threshold %in% Thresholds), aes(x=log(mmHgH), fill=as.factor(Mort6mo))) +
+  geom_histogram() +
+  scale_fill_manual(values=c("deepskyblue", "orange3"), name="6 month mortality") +
+  facet_wrap(~Threshold, scales = "free") +
+  theme_minimal() + 
+  theme(legend.title = element_text(face="italic"), text = element_text(size = 14), axis.text = element_text(size=14))
+
+dev.new()
+print(distr)
+
+ FAV.distr <- ggplot(data=PTD$tot$PTD.long %>% filter(Threshold %in% Thresholds), aes(x=mmHgH.perday, fill=as.factor(Outcome))) +
+  geom_histogram() +
+  scale_fill_manual(values=c("deepskyblue", "orange3"), name="Outcome") +
+  facet_wrap(~Threshold, scales = "free") +
+  theme_minimal() + 
+  theme(legend.title = element_text(face="italic"), text = element_text(size = 14), axis.text = element_text(size=14))
+
+# dev.new()
+# print(FAV.distr)
+
+# dev.new()
+# print(ggplot(PTD$tot$PTD.long %>% filter(Threshold %in% Thresholds), aes(x=mmHgH)) + geom_histogram() + facet_wrap(~Threshold, scales="free_x"))
+
+
+
+
+
+# PTD stratified per GOSE:
+jet.colors <- colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F", "yellow", "#FF7F00", "red"))#, "#7F0000"))
+# library(RColorBrewer)
+
+dev.new()
+print(ggplot(data=PTD$tot$PTD.GOSE, aes(Threshold, mean.mmHgH, group=GOSE6mo, color=GOSE6mo)) + geom_line(size=2) +
+        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=14), axis.text=element_text(size = 14)) +
+        #xlim(0,30) +
+        ggtitle("Mean PTD stratified by 6-month outcome") + xlab("ICP threshold (mmHg)") + ylab("mean mmHg"))
+
+
+
+# Fig Mean PTD, favorable/unfavorable outcome:
+
+PTD.fav.plot <- ggplot(PTD$tot$PTD.FAV, aes(x=Threshold, y=mean.mmHgH, color=Outcome)) + geom_line(size=2) +
+        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"), text = element_text(size = 20), axis.text = element_text(size=20)) +
+        ggtitle("A. Favorable vs Unfavorable outcome") + xlab("ICP threshold") + ylab("PTD (Mean mmHg)")
+
+PTD.mort.plot <- ggplot(PTD$tot$PTD.Mort, aes(x=Threshold, y=mean.mmHgH, color=Mort6mo)) + geom_line(size=2) +
+        scale_color_manual(values=c("deepskyblue", "orange3"), name="6 month mortality") +
+        #labs(color = "6 month mortality") +
+        theme_minimal() +
+        theme(legend.title = element_text(face="italic"), text = element_text(size = 20), axis.text = element_text(size=14)) +
+        ggtitle("B. 6 month mortality") + xlab("ICP threshold") + ylab("PTD (Mean mmHg)")
+
+dev.new()
+grid.arrange(PTD.fav.plot, PTD.mort.plot, ncol=2)
+
+# Fig active/passive (not normalized):
+Active.Fav.plot <- ggplot(PTD$active$PTD.FAV, aes(x=Threshold, y=mean.mmHgH, color=Outcome)) +
+  geom_line(size=2) +
+  scale_color_manual(values=c("deepskyblue", "orange3"), name="Outcome") +
+  theme_minimal() +
+  theme(legend.title = element_text(face="italic"), text=element_text(size=14), axis.text = element_text(size = 14)) +
+#  ylim(0,25) +
+  ggtitle("A. Intact autoregulation") + xlab("ICP threshold") + ylab("PTD (Mean mmHg*h/time \n of intact autoregulation)")
+
+Passive.Fav.plot <- ggplot(PTD$passive$PTD.FAV, aes(x=Threshold, y=mean.mmHgH, color=Outcome)) +
+  geom_line(size=2) +
+  scale_color_manual(values=c("deepskyblue", "orange3"), name="Outcome") +
+  theme_minimal() +
+  theme(legend.title = element_text(face="italic"), text=element_text(size=14), axis.text = element_text(size = 14)) +
+ # ylim(0,25) +
+  ggtitle("B. Impaired autoregulation") + xlab("ICP threshold") + ylab("PTD (Mean mmHg*h/time \n of impaired autoregulation)")
+
+Active.Mort.plot <- ggplot(PTD$active$PTD.Mort, aes(x=Threshold, y=mean.mmHgH, color=Mort6mo)) +
+  geom_line(size=2) +
+  scale_color_manual(values=c("deepskyblue", "orange3"), name="6 month mortality") +
+  theme_minimal() +
+  theme(legend.title = element_text(face="italic"), text=element_text(size=14), axis.text = element_text(size = 14)) +
+  #ylim(0,25) +
+  ggtitle("C. Intact autoregulation") + xlab("ICP threshold") + ylab("PTD (Mean mmHg*h/time \n of intact autoregulation)")
+
+Passive.Mort.plot <- ggplot(PTD$passive$PTD.Mort, aes(x=Threshold, y=mean.mmHgH, color=Mort6mo)) +
+  geom_line(size=2) +
+  scale_color_manual(values=c("deepskyblue", "orange3"), name="6 month mortality") +
+  theme_minimal() +
+  theme(legend.title = element_text(face="italic"), text=element_text(size=14), axis.text = element_text(size = 14)) +
+  #ylim(0,25) +
+  ggtitle("D. Impaired autoregulation") + xlab("ICP threshold") + ylab("PTD (Mean mmHg*h/time \n of impaired autoregulation)")
+
+
+library(grid)
+
+dev.new()
+grid.arrange(arrangeGrob(Active.Fav.plot, Passive.Fav.plot, top = textGrob("Favorable vs Unfavorable outcome", gp=gpar(fontsize=20)), ncol=2),
+             arrangeGrob(Active.Mort.plot, Passive.Mort.plot, top=textGrob("6-month mortality", gp=gpar(fontsize=20)), ncol=2))
+
+
+# Tables: ----------------------------------------------------------
+
+
+# Kolmogorov-Smirnov: ------------------------
+KS.Table.fcn <- function(df, decimals) # %>% select(gupi:mmHgH.perday, Outcome) %>% rename(Outcome = )
+{
+  x <- df %>% filter(Threshold == Thresholds[1]) %>% group_by(Outcome) %>% summarise(Mean = round(mean(dose),1),
+                                                                                     SD=round(sd(dose),1))
+                                                                                     
+  x$p <- round(ks.test(x=filter(df, Threshold==Thresholds[1], Outcome==x$Outcome[1])$dose, y=filter(df, Threshold==Thresholds[1], Outcome==x$Outcome[2])$dose)$p.value,2)
+  
+  #x$p <- round(ks.test(data=df %>% filter(Threshold == Thresholds[1]), dose~Outcome)$p.value, 3)
+  x$PTD <- Thresholds[1]
+  
+  for(i in 2:length(Thresholds))
+  {
+    temp <- df %>% filter(Threshold == Thresholds[i]) %>% group_by(Outcome) %>% summarise(Mean = round(mean(dose),1),
+                                                                                          SD=round(sd(dose),1))
+    temp$p <- round(ks.test(x=filter(df, Threshold==Thresholds[i], Outcome==x$Outcome[1])$dose, y=filter(df, Threshold==Thresholds[i], Outcome==x$Outcome[2])$dose)$p.value, decimals)
+    temp$PTD <- Thresholds[i]
+    x <- rbind(x, temp)
+  }
+  
+  x <- x %>% mutate(Mean = paste0(Mean, " (±", SD, ")")) %>% select(-SD)
+  
+  x <- spread(x, Outcome, Mean)
+  x <- x %>% arrange(PTD) %>% select(-p, PTD, everything(), p)
+  
+  return(x)
+}
+
+# Table 2: PTD over certain thresholds ---------------------------
+Tot.Fav <- KS.Table.fcn(df = PTD$tot$PTD.long %>% rename(dose=mmHgH), decimals=3)
+Tot.Fav <- Tot.Fav %>% rename("Favorable outcome" = FAV, "Unfavorable outcome" = UNFAV)
+
+Tot.Mort <- KS.Table.fcn(df=PTD$tot$PTD.long %>% select(Threshold, mmHgH, Mort6mo) %>% rename(Outcome=Mort6mo, dose=mmHgH), decimals=3)
+Tot.Mort <- Tot.Mort %>% rename("Dead at 6 mo" = Yes, "Alive at 6 mo" = No)
+
+View(Tot.Fav)
+View(Tot.Mort)
+
+# Un-normalized PTD:s:
+# Table 3: PTD by autoregulation status, Favorable/Unfavorable: ---------
+Active.Fav <- KS.Table.fcn(df=PTD$active$PTD.long %>% select(gupi:Threshold, Outcome, mmHgH) %>% rename(dose=mmHgH), decimals=3)
+Active.Fav <- Active.Fav %>% rename("Favorable Outcome" = FAV, "Unfavorable outcome" = UNFAV)
+
+Passive.Fav <- KS.Table.fcn(df=PTD$passive$PTD.long %>% select(gupi:Threshold, Outcome, mmHgH) %>% rename(dose=mmHgH), decimals=3)
+Passive.Fav <- Passive.Fav %>% rename("Favorable Outcome" = FAV, "Unfavorable outcome" = UNFAV)
+
+View(Active.Fav)
+View(Passive.Fav)
+
+# Table 4:
+Active.Mort <- KS.Table.fcn(df=PTD$active$PTD.long %>% select(gupi:Threshold, Mort6mo, mmHgH) %>% rename(dose=mmHgH, Outcome = Mort6mo), decimals=3)
+Active.Mort <- Active.Mort %>% rename("Dead at 6 mo" = Yes, "Alive at 6 mo" = No)
+
+Passive.Mort <- KS.Table.fcn(df=PTD$passive$PTD.long %>% select(gupi:Threshold, Mort6mo, mmHgH) %>% rename(dose=mmHgH, Outcome = Mort6mo), decimals=3)
+Passive.Mort <- Passive.Mort %>% rename("Dead at 6 mo" = Yes, "Alive at 6 mo" = No)
+
+View(Active.Mort)
+View(Passive.Mort)
+
+# ----- THE END ------------------------------------------------
+
+
+
+# Time over thresholds: -------------------------------------------
+View(ms %>% filter(gupi == ms$gupi[1]))
+
+Times <- function(i, dta){dta %>%  group_by(gupi) %>% filter(icp >= i) %>% summarise(minutes = n())}
+TimeOverThresholds <- list()
+
+TimeOverThresholds$tot <- lapply(0:40, Times, ms)
+TimeOverThresholds$passive <- lapply(0:40, Times, ms %>% filter(PRx >= 0.3))
+TimeOverThresholds$active <- lapply(0:40, Times, ms %>% filter(PRx < 0.3))
+
+for(j in 1:length(TimeOverThresholds))
+{
+  for (i in 1:length(TimeOverThresholds[[j]]))
+  {
+    TimeOverThresholds[[j]][[i]] <- TimeOverThresholds[[j]][[i]] %>% mutate(Threshold = i-1)
+  }
+}
+
+View(TimeOverThresholds$tot)
+
+for(i in 1:length(TimeOverThresholds))
+{
+  TimeOverThresholds[[i]] <- bind_rows(TimeOverThresholds[[i]]) %>% mutate(hours=round(minutes/60, 2)) %>% select(-minutes)
+  TimeOverThresholds[[i]] <- left_join(TimeOverThresholds[[i]], select(Baseline, gupi, GOSE6mo, GOS, GOSE6mo.derived, Mort90d, Mort6mo, Outcome), by="gupi")
+}
+
+
+View(TimeOverThresholds$tot)
+
+TOT.summaries <- list()
+Temp <- list()
+
+for(i in 1:length(TimeOverThresholds))
+{
+  Temp$GOSE <- TimeOverThresholds[[i]] %>% group_by(GOSE6mo, Threshold) %>% summarise(mean.time = sum(hours)/length(unique(ms$gupi)))
+Temp$FAV <- TimeOverThresholds[[i]] %>% group_by(Outcome, Threshold) %>% summarise(mean.time = sum(hours)/length(unique(ms$gupi)))
+ Temp$Mort <- TimeOverThresholds[[i]] %>% group_by(Mort6mo, Threshold) %>% summarise(mean.time = sum(hours)/length(unique(ms$gupi)))
+  TOT.summaries[[i]] <- Temp
+}
+names(TOT.summaries) <- names(TimeOverThresholds)
+
+
+
+dev.new()
+print(ggplot(data=TOT.summaries$tot$GOSE, aes(Threshold, mean.time, group=GOSE6mo, color=GOSE6mo)) + geom_line(size=2) +
+        scale_color_manual(values=rev(jet.colors(8)), name="6 month GOS-E") +
+        theme_minimal() +
+        theme(legend.title = element_text(face="italic"), text=element_text(size=14), axis.text=element_text(size = 14)) +
+        #xlim(20,30) + ylim(0,6) +
+        ggtitle("Time over thresholds") + xlab("ICP threshold (mmHg)") + ylab("mean time (hours)"))
+
+dev.new()
+print(ggplot(data=TOT.summaries$tot$FAV, aes(Threshold, mean.time, group=Outcome, color=Outcome)) + geom_line(size=2) +
+        scale_color_manual(values=rev(jet.colors(8)), name="Favorable/Unfavorable outcome") +
+        theme_minimal() +
+        theme(legend.title = element_text(face="italic"), text=element_text(size=14), axis.text=element_text(size = 14)) +
+        xlim(0,30) +
+        ggtitle("Time over thresholds") + xlab("ICP threshold (mmHg)") + ylab("mean time (hours)"))
+
+dev.new()
+print(ggplot(data=TOT.summaries$tot$Mort, aes(Threshold, mean.time, group=Mort6mo, color=Mort6mo)) + geom_line(size=2) +
+        scale_color_manual(values=rev(jet.colors(8)), name="6 month mortality") +
+        theme_minimal() +
+        theme(legend.title = element_text(face="italic"), text=element_text(size=14), axis.text=element_text(size = 14)) +
+        #xlim(20,30) + ylim(0,6) +
+        ggtitle("Time over thresholds") + xlab("ICP threshold (mmHg)") + ylab("mean time (hours)"))
+
+
+
+# Data to Anders: -------------------------------------------------
+
+AUC.tot.temp <- left_join(AUC.tot %>% select(-contains("perday")), select(Baseline, gupi, at, pt, monitor.time))
+AUC.tot.temp$ID <- 1:nrow(AUC.tot.temp)
+AUC.tot.temp <- AUC.tot.temp %>% select(ID, everything(), -gupi)
+
+write.csv(AUC.tot.temp, file=paste0(f,"200318/ICP.AUC.csv"), row.names=FALSE)
+
+View(AUC.tot.temp)
+
+AUC.active.norm.temp <- AUC.active.norm %>% select(-contains("perhour"))
+AUC.active.norm.temp$ID <- 1:nrow(AUC.active.norm.temp)
+AUC.active.norm.temp <- AUC.active.norm.temp %>% select(ID, everything(), -gupi)
+write.csv(AUC.active.norm.temp, file=paste0(f,"200318/ICP.AUC.active.csv"), row.names=FALSE)
+
+AUC.passive.norm.temp <- AUC.passive.norm %>% select(-contains("perhour"))
+AUC.passive.norm.temp$ID <- 1:nrow(AUC.passive.norm.temp)
+AUC.passive.norm.temp <- AUC.passive.norm.temp %>% select(ID, everything(), -gupi)
+
+View(AUC.passive.norm.temp)
+write.csv(AUC.passive.norm.temp, file=paste0(f,"200318/ICP.AUC.passive.csv"), row.names = FALSE)
+
+