Upload New File

ICP_burden_multivar_201006.R

+rm(list=ls())
+
+library(dplyr)
+library(tidyr)
+library(ggplot2)
+library(gridExtra)
+
+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_preparation_200225.csv
+
+
+Baseline  <- Baseline %>% mutate(Mort6mo = ifelse(GOSE6mo == 1, 1, 0),
+                                 Outcome = as.factor(ifelse(GOSE6mo %in% c(1,2,3,4), "UNFAV", "FAV")),
+                                 icp.mean = ifelse(is.na(icp.mean), icp.evd.mean, icp.mean),
+                                 icp.median = ifelse(is.na(icp.median), icp.evd.median, icp.median),
+                                 cpp.mean = ifelse(is.na(cpp.mean), cpp.evd.mean, cpp.mean),
+                                 cpp.median = ifelse(is.na(cpp.median), cpp.evd.median, cpp.median))
+
+load(file=paste0(f,"200625/PTD.RData"))
+
+# TIL.max added as predictor: ---------------------------------------------------------------------------------------
+# Logistic regression model for time in red zone: ----
+
+# Univariate regression:
+uni.rz.outcome <- glm(Outcome~rz, family=binomial(), data=Baseline)
+exp(uni.rz.outcome$coefficients) #OR
+exp(confint(uni.rz.outcome)) # confidence interval 95%
+summary(uni.rz.outcome)$coefficients[2,4] # p value
+
+uni.rz.mort <- glm(Mort6mo~rz, family=binomial(), data=Baseline)
+exp(uni.rz.mort$coefficients)
+exp(confint(uni.rz.mort))
+summary(uni.rz.mort)$coefficients[2,4]
+
+
+# 6 month mortality / unfav/fav outcome:
+mort.model <- glm(data=Baseline, Mort6mo~rz + Age + GCSMotorBaselineDerived + PupilsBaseline + TIL.max, family=binomial())
+fav.model <- glm(Outcome~rz + Age + GCSMotorBaselineDerived + PupilsBaseline + TIL.max, family="binomial", data=Baseline)
+
+summary(mort.model)
+summary(fav.model)
+
+
+OR.mort <- round(as.data.frame(exp(cbind(OR=coef(mort.model), confint(mort.model)))),2)
+OR.mort$p <- round(summary(mort.model)$coefficients[,4],3)
+OR.mort <- OR.mort %>% mutate(Variable=rownames(OR.mort), "OR, 95% CI" = paste0(OR, " (",`2.5 %`, "-", `97.5 %`,")")) %>% select(Variable, "OR, 95% CI", p)
+
+OR.fav <- round(as.data.frame(exp(cbind(OR=coef(fav.model), confint(fav.model)))),2)
+OR.fav$p <- round(summary(fav.model)$coefficients[,4],3)
+OR.fav <- OR.fav %>% mutate(Variable=rownames(OR.fav), "OR, 95% CI" = paste0(OR, " (",`2.5 %`, "-", `97.5 %`,")")) %>% select(Variable, "OR, 95% CI", p)
+
+View(OR.mort)
+View(OR.fav)
+
+
+# Time in red zone above +2sd:
+mort.model <- glm(data=Baseline, Mort6mo~rz.ub + Age + GCSMotorBaselineDerived + PupilsBaseline + TIL.max, family=binomial())
+fav.model <- glm(Outcome~rz.ub + Age + GCSMotorBaselineDerived + PupilsBaseline + TIL.max, family="binomial", data=Baseline)
+
+OR.mort.ub <- round(as.data.frame(exp(cbind(OR=coef(mort.model), confint(mort.model)))),2)
+OR.mort.ub$p <- round(summary(mort.model)$coefficients[,4],3)
+OR.mort.ub <- OR.mort.ub %>% mutate(Variable=rownames(OR.mort.ub), "OR, 95% CI" = paste0(OR, " (",`2.5 %`, "-", `97.5 %`,")")) %>% select(Variable, "OR, 95% CI", p)
+
+OR.fav.ub <- round(as.data.frame(exp(cbind(OR=coef(fav.model), confint(fav.model)))),2)
+OR.fav.ub$p <- round(summary(fav.model)$coefficients[,4],3)
+OR.fav.ub <- OR.fav.ub %>% mutate(Variable=rownames(OR.fav.ub), "OR, 95% CI" = paste0(OR, " (",`2.5 %`, "-", `97.5 %`,")")) %>% select(Variable, "OR, 95% CI", p)
+
+View(OR.mort.ub)
+View(OR.fav.ub)
+
+
+# Multivariate logistic regression model: ----------------------------------
+
+# Add to baseline, PTD active/passive/tot at thresholds
+Thresholds <- c(0,10,15,20,25,30)
+# View(PTD$tot$PTD.long)
+PTD.tot <- spread(PTD$tot$PTD.long %>% select(gupi, Threshold, mmHgH) %>% filter(Threshold %in% Thresholds), Threshold, mmHgH)
+colnames(PTD.tot)[2:7] <- paste0("PTD.tot.",colnames(PTD.tot[2:7]))
+
+PTD.active <- spread(PTD$active$PTD.long %>% select(gupi, Threshold, mmHgH) %>% filter(Threshold %in% Thresholds), Threshold, mmHgH)
+colnames(PTD.active)[2:7] <- paste0("PTD.active.",colnames(PTD.active[2:7]))
+PTD.passive <- spread(PTD$passive$PTD.long %>% select(gupi, Threshold, mmHgH) %>% filter(Threshold %in% Thresholds), Threshold, mmHgH)
+colnames(PTD.passive)[2:7] <- paste0("PTD.passive.",colnames(PTD.passive[2:7]))
+
+Baseline <- left_join(Baseline, PTD.tot, by="gupi")
+Baseline <- left_join(Baseline, PTD.active, by="gupi")
+Baseline <- left_join(Baseline, PTD.passive, by="gupi")
+
+
+# OR: ---------------
+
+OR.fcn <- function(x)
+{
+  y <- round(as.data.frame(exp(cbind(OR=coef(x), conf.int=confint(x)))), 5)
+  y <- y %>% mutate(Variable = rownames(y), p=round(summary(x)$coefficients[,4],3))
+  y <- y %>% mutate("OR, 95% CI" = paste0(OR, ", (", `2.5 %`, "-", `97.5 %`, ")")) %>% select(Variable, "OR, 95% CI", p)
+  return(y)
+}
+
+OR.list <- list()
+exp(summary(glm(Mort6mo~PTD.tot.0, data=Baseline, family=binomial()))$coef)
+mean((Baseline %>% filter(Mort6mo==0))$PTD.tot.0, na.rm=TRUE)
+
+# Total PTD:
+OR.list$Mort6mo$PTD.tot$'0' <- OR.fcn(glm(Mort6mo~PTD.tot.0 + Age + GCSMotorBaselineDerived + PupilsBaseline + TIL.max, family="binomial", data=Baseline))
+OR.list$Mort6mo$PTD.tot$'10' <- OR.fcn(glm(Mort6mo~PTD.tot.10 + Age + GCSMotorBaselineDerived + PupilsBaseline + TIL.max, family="binomial", data=Baseline))
+OR.list$Mort6mo$PTD.tot$'15' <- OR.fcn(glm(Mort6mo~PTD.tot.15 + Age + GCSMotorBaselineDerived + PupilsBaseline + TIL.max, family="binomial", data=Baseline))
+OR.list$Mort6mo$PTD.tot$'20' <- OR.fcn(glm(Mort6mo~PTD.tot.20 + Age + GCSMotorBaselineDerived + PupilsBaseline + TIL.max, family="binomial", data=Baseline))
+OR.list$Mort6mo$PTD.tot$'25' <- OR.fcn(glm(Mort6mo~PTD.tot.25 + Age + GCSMotorBaselineDerived + PupilsBaseline + TIL.max, family="binomial", data=Baseline))
+OR.list$Mort6mo$PTD.tot$'30' <- OR.fcn(glm(Mort6mo~PTD.tot.30 + Age + GCSMotorBaselineDerived + PupilsBaseline + TIL.max, family="binomial", data=Baseline))
+
+# Active PTD:
+OR.list$Mort6mo$PTD.active$'0' <- OR.fcn(glm(Mort6mo~PTD.active.0 + Age + GCSMotorBaselineDerived + PupilsBaseline + TIL.max, family="binomial", data=Baseline))
+OR.list$Mort6mo$PTD.active$'10' <- OR.fcn(glm(Mort6mo~PTD.active.10 + Age + GCSMotorBaselineDerived + PupilsBaseline + TIL.max, family="binomial", data=Baseline))
+OR.list$Mort6mo$PTD.active$'15' <- OR.fcn(glm(Mort6mo~PTD.active.15 + Age + GCSMotorBaselineDerived + PupilsBaseline + TIL.max, family="binomial", data=Baseline))
+OR.list$Mort6mo$PTD.active$'20' <- OR.fcn(glm(Mort6mo~PTD.active.20 + Age + GCSMotorBaselineDerived + PupilsBaseline + TIL.max, family="binomial", data=Baseline))
+OR.list$Mort6mo$PTD.active$'25' <- OR.fcn(glm(Mort6mo~PTD.active.25 + Age + GCSMotorBaselineDerived + PupilsBaseline + TIL.max, family="binomial", data=Baseline))
+OR.list$Mort6mo$PTD.active$'30' <- OR.fcn(glm(Mort6mo~PTD.active.30 + Age + GCSMotorBaselineDerived + PupilsBaseline + TIL.max, family="binomial", data=Baseline))
+
+# Passive PTD:
+OR.list$Mort6mo$PTD.passive$'0' <- OR.fcn(glm(Mort6mo~PTD.passive.0 + Age + GCSMotorBaselineDerived + PupilsBaseline + TIL.max, family="binomial", data=Baseline))
+OR.list$Mort6mo$PTD.passive$'10' <- OR.fcn(glm(Mort6mo~PTD.passive.10 + Age + GCSMotorBaselineDerived + PupilsBaseline + TIL.max, family="binomial", data=Baseline))
+OR.list$Mort6mo$PTD.passive$'15' <- OR.fcn(glm(Mort6mo~PTD.passive.15 + Age + GCSMotorBaselineDerived + PupilsBaseline + TIL.max, family="binomial", data=Baseline))
+OR.list$Mort6mo$PTD.passive$'20' <- OR.fcn(glm(Mort6mo~PTD.passive.20 + Age + GCSMotorBaselineDerived + PupilsBaseline + TIL.max, family="binomial", data=Baseline))
+OR.list$Mort6mo$PTD.passive$'25' <- OR.fcn(glm(Mort6mo~PTD.passive.25 + Age + GCSMotorBaselineDerived + PupilsBaseline + TIL.max, family="binomial", data=Baseline))
+OR.list$Mort6mo$PTD.passive$'30' <- OR.fcn(glm(Mort6mo~PTD.passive.30 + Age + GCSMotorBaselineDerived + PupilsBaseline + TIL.max, family="binomial", data=Baseline))
+
+
+# Total PTD:
+OR.list$Outcome$PTD.tot$'0' <- OR.fcn(glm(Outcome~PTD.tot.0 + Age + GCSMotorBaselineDerived + PupilsBaseline + TIL.max, family="binomial", data=Baseline))
+OR.list$Outcome$PTD.tot$'10' <- OR.fcn(glm(Outcome~PTD.tot.10 + Age + GCSMotorBaselineDerived + PupilsBaseline + TIL.max, family="binomial", data=Baseline))
+OR.list$Outcome$PTD.tot$'15' <- OR.fcn(glm(Outcome~PTD.tot.15 + Age + GCSMotorBaselineDerived + PupilsBaseline + TIL.max, family="binomial", data=Baseline))
+OR.list$Outcome$PTD.tot$'20' <- OR.fcn(glm(Outcome~PTD.tot.20 + Age + GCSMotorBaselineDerived + PupilsBaseline + TIL.max, family="binomial", data=Baseline))
+OR.list$Outcome$PTD.tot$'25' <- OR.fcn(glm(Outcome~PTD.tot.25 + Age + GCSMotorBaselineDerived + PupilsBaseline + TIL.max, family="binomial", data=Baseline))
+OR.list$Outcome$PTD.tot$'30' <- OR.fcn(glm(Outcome~PTD.tot.30 + Age + GCSMotorBaselineDerived + PupilsBaseline + TIL.max, family="binomial", data=Baseline))
+
+# Active PTD:
+OR.list$Outcome$PTD.active$'0' <- OR.fcn(glm(Outcome~PTD.active.0 + Age + GCSMotorBaselineDerived + PupilsBaseline + TIL.max, family="binomial", data=Baseline))
+OR.list$Outcome$PTD.active$'10' <- OR.fcn(glm(Outcome~PTD.active.10 + Age + GCSMotorBaselineDerived + PupilsBaseline + TIL.max, family="binomial", data=Baseline))
+OR.list$Outcome$PTD.active$'15' <- OR.fcn(glm(Outcome~PTD.active.15 + Age + GCSMotorBaselineDerived + PupilsBaseline + TIL.max, family="binomial", data=Baseline))
+OR.list$Outcome$PTD.active$'20' <- OR.fcn(glm(Outcome~PTD.active.20 + Age + GCSMotorBaselineDerived + PupilsBaseline + TIL.max, family="binomial", data=Baseline))
+OR.list$Outcome$PTD.active$'25' <- OR.fcn(glm(Outcome~PTD.active.25 + Age + GCSMotorBaselineDerived + PupilsBaseline + TIL.max, family="binomial", data=Baseline))
+OR.list$Outcome$PTD.active$'30' <- OR.fcn(glm(Outcome~PTD.active.30 + Age + GCSMotorBaselineDerived + PupilsBaseline + TIL.max, family="binomial", data=Baseline))
+
+# Passive PTD:
+OR.list$Outcome$PTD.passive$'0' <- OR.fcn(glm(Outcome~PTD.passive.0 + Age + GCSMotorBaselineDerived + PupilsBaseline + TIL.max, family="binomial", data=Baseline))
+OR.list$Outcome$PTD.passive$'10' <- OR.fcn(glm(Outcome~PTD.passive.10 + Age + GCSMotorBaselineDerived + PupilsBaseline + TIL.max, family="binomial", data=Baseline))
+OR.list$Outcome$PTD.passive$'15' <- OR.fcn(glm(Outcome~PTD.passive.15 + Age + GCSMotorBaselineDerived + PupilsBaseline + TIL.max, family="binomial", data=Baseline))
+OR.list$Outcome$PTD.passive$'20' <- OR.fcn(glm(Outcome~PTD.passive.20 + Age + GCSMotorBaselineDerived + PupilsBaseline + TIL.max, family="binomial", data=Baseline))
+OR.list$Outcome$PTD.passive$'25' <- OR.fcn(glm(Outcome~PTD.passive.25 + Age + GCSMotorBaselineDerived + PupilsBaseline + TIL.max, family="binomial", data=Baseline))
+OR.list$Outcome$PTD.passive$'30' <- OR.fcn(glm(Outcome~PTD.passive.30 + Age + GCSMotorBaselineDerived + PupilsBaseline + TIL.max, family="binomial", data=Baseline))
+
+
+# Table C, supplement:
+View(OR.list$Outcome$PTD.tot$`20`)
+View(OR.list$Mort6mo$PTD.tot$`20`)
+
+# Table D, supplement:
+View(OR.list$Outcome$PTD.active$`20`)
+View(OR.list$Mort6mo$PTD.active$`20`)
+
+# Table E, supplement:
+View(OR.list$Outcome$PTD.passive$`20`)
+View(OR.list$Mort6mo$PTD.passive$`20`)
+
+