#  Important packages to load in
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(gtsummary)
library(survival)
library(labelled)
library(tidyverse)
## ── Attaching packages ──────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.3     ✓ purrr   0.3.4
## ✓ tibble  3.0.3     ✓ stringr 1.4.0
## ✓ tidyr   1.1.1     ✓ forcats 0.5.0
## ✓ readr   1.3.1
## ── Conflicts ─────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(broom)
library(ggplot2)
library(survminer)
## Loading required package: ggpubr
library(VIM)
## Loading required package: colorspace
## Loading required package: grid
## VIM is ready to use.
## Suggestions and bug-reports can be submitted at: https://github.com/statistikat/VIM/issues
## 
## Attaching package: 'VIM'
## The following object is masked from 'package:datasets':
## 
##     sleep
library(tibble)
library(ggplot2)
library(OddsPlotty)
library(e1071)
library(ggthemes)
library(mlbench)
library(caret)
## Loading required package: lattice
## 
## Attaching package: 'caret'
## The following object is masked from 'package:purrr':
## 
##     lift
## The following object is masked from 'package:survival':
## 
##     cluster
load("TextbookOutcomes.RData")
TextbookOutcomes <- 
  TextbookOutcomes %>%
  mutate(TextbookOutcomes, Age = ifelse(Age >= 85, "85+", 
                                               ifelse(Age >= 75, "75-84", 
                                                      ifelse(Age >= 65, "65-74", 
                                                             ifelse(Age >=0, "0-65", NA)))))
### Recoding mortality


TextbookOutcomes$DOpertoD[TextbookOutcomes$DOpertoD == -99] <- "No"

TextbookOutcomes$DOpertoD[TextbookOutcomes$DOpertoD == 0] <- "Yes"
TextbookOutcomes$DOpertoD[TextbookOutcomes$DOpertoD == 1] <- "Yes"
TextbookOutcomes$DOpertoD[TextbookOutcomes$DOpertoD == 2] <- "Yes"
TextbookOutcomes$DOpertoD[TextbookOutcomes$DOpertoD == 3] <- "Yes"
TextbookOutcomes$DOpertoD[TextbookOutcomes$DOpertoD == 4] <- "Yes"
TextbookOutcomes$DOpertoD[TextbookOutcomes$DOpertoD == 5] <- "Yes"
TextbookOutcomes$DOpertoD[TextbookOutcomes$DOpertoD == 6] <- "Yes"
TextbookOutcomes$DOpertoD[TextbookOutcomes$DOpertoD == 7] <- "Yes"
TextbookOutcomes$DOpertoD[TextbookOutcomes$DOpertoD == 8] <- "Yes"
TextbookOutcomes$DOpertoD[TextbookOutcomes$DOpertoD == 9] <- "Yes"
TextbookOutcomes$DOpertoD[TextbookOutcomes$DOpertoD == 10] <- "Yes"
TextbookOutcomes$DOpertoD[TextbookOutcomes$DOpertoD == 11] <- "Yes"
TextbookOutcomes$DOpertoD[TextbookOutcomes$DOpertoD == 12] <- "Yes"
TextbookOutcomes$DOpertoD[TextbookOutcomes$DOpertoD == 13] <- "Yes"
TextbookOutcomes$DOpertoD[TextbookOutcomes$DOpertoD == 14] <- "Yes"
TextbookOutcomes$DOpertoD[TextbookOutcomes$DOpertoD == 15] <- "Yes"
TextbookOutcomes$DOpertoD[TextbookOutcomes$DOpertoD == 16] <- "Yes"
TextbookOutcomes$DOpertoD[TextbookOutcomes$DOpertoD == 17] <- "Yes"
TextbookOutcomes$DOpertoD[TextbookOutcomes$DOpertoD == 18] <- "Yes"
TextbookOutcomes$DOpertoD[TextbookOutcomes$DOpertoD == 19] <- "Yes"
TextbookOutcomes$DOpertoD[TextbookOutcomes$DOpertoD == 20] <- "Yes"
TextbookOutcomes$DOpertoD[TextbookOutcomes$DOpertoD == 21] <- "Yes"
TextbookOutcomes$DOpertoD[TextbookOutcomes$DOpertoD == 22] <- "Yes"
TextbookOutcomes$DOpertoD[TextbookOutcomes$DOpertoD == 23] <- "Yes"
TextbookOutcomes$DOpertoD[TextbookOutcomes$DOpertoD == 24] <- "Yes"
TextbookOutcomes$DOpertoD[TextbookOutcomes$DOpertoD == 25] <- "Yes"
TextbookOutcomes$DOpertoD[TextbookOutcomes$DOpertoD == 26] <- "Yes"
TextbookOutcomes$DOpertoD[TextbookOutcomes$DOpertoD == 27] <- "Yes"
TextbookOutcomes$DOpertoD[TextbookOutcomes$DOpertoD == 28] <- "Yes"
TextbookOutcomes$DOpertoD[TextbookOutcomes$DOpertoD == 29] <- "Yes"
TextbookOutcomes$DOpertoD[TextbookOutcomes$DOpertoD == 30] <- "Yes"
TextbookOutcomes$PRNCPTX[TextbookOutcomes$PRNCPTX == "EXCISION AMPULLA VATER"] <- NA
TextbookOutcomes$PRNCPTX[TextbookOutcomes$PRNCPTX == "EXCISION LESION PANCREAS"] <- NA
TextbookOutcomes$PRNCPTX[TextbookOutcomes$PRNCPTX == "UNLISTED PROCEDURE PANCREAS"] <- NA
TextbookOutcomes$PRNCPTX[TextbookOutcomes$PRNCPTX == "PNCRTECT WHIPPLE W/O PANCREATOJEJUNOSTOMY"] <- "Whipple"
TextbookOutcomes$PRNCPTX[TextbookOutcomes$PRNCPTX == "PNCRTECT W/PANCREATOJEJUNOSTOMY"] <- "Whipple"
TextbookOutcomes$PRNCPTX[TextbookOutcomes$PRNCPTX == "PNCRTECT DSTL STOT W/PNCRTCOJEJUNOSTOMY"] <- "Distal Pancreatectomy"
TextbookOutcomes$PRNCPTX[TextbookOutcomes$PRNCPTX == "PNCRTECT DSTL STOT W/O PNCRTCOJEJUNOSTOMY"] <- "Distal Pancreatectomy"
TextbookOutcomes$PRNCPTX[TextbookOutcomes$PRNCPTX == "PNCRTECT PROX STOT W/PANCREATOJEJUNOSTOMY"] <- "Whipple"
TextbookOutcomes$PRNCPTX[TextbookOutcomes$PRNCPTX == "PNCRTECT PROX STOT W/O PANCREATOJEJUNOSTOMY"] <- "Whipple"
TextbookOutcomes$PRNCPTX[TextbookOutcomes$PRNCPTX == "PANCREATECTOMY TOTAL"] <- "Total Pancreatectomy"
TextbookOutcomes$PRNCPTX[TextbookOutcomes$PRNCPTX == "PNCRTECT DSTL NR-TOT W/PRSRV DUO CHLD-TYP PX"] <- "Total Pancreatectomy"


TextbookOutcomes$PRNCPTX[TextbookOutcomes$PRNCPTX == "HEPATECTOMY RESCJ PARTIAL LOBECTOMY"] <- "Partial Lobectomy"
TextbookOutcomes$PRNCPTX[TextbookOutcomes$PRNCPTX == "HEPATECTOMY RESCJ TOTAL LEFT LOBECTOMY"] <- "Lobectomy"
TextbookOutcomes$PRNCPTX[TextbookOutcomes$PRNCPTX == "HEPATECTOMY RESCJ TOTAL RIGHT LOBECTOMY"] <- "Lobectomy"
TextbookOutcomes$PRNCPTX[TextbookOutcomes$PRNCPTX == "HEPATECTOMY RESCJ TRISEGMENTECTOMY"] <- "Trisegmentectomy"


table(TextbookOutcomes$PRNCPTX)
## 
## Distal Pancreatectomy             Lobectomy     Partial Lobectomy 
##                  8292                  4274                 11817 
##  Total Pancreatectomy      Trisegmentectomy               Whipple 
##                   706                  1366                 16464
TextbookOutcomes <-
  TextbookOutcomes %>%
  mutate(TextbookOutcomes, PatientGroup = ifelse(PRNCPTX == "Partial Lobectomy", "Liver",
                                                 ifelse(PRNCPTX == "Lobectomy", "Liver",
                                                        ifelse(PRNCPTX == "Trisegmentectomy", "Liver",
                                                               ifelse(PRNCPTX == "Whipple", "Pancreas",
                                                                             ifelse(PRNCPTX == "Distal Pancreatectomy", "Pancreas",
                                                                                    ifelse(PRNCPTX == "Total Pancreatectomy", "Pancreas", NA)))))))
### Recoding readmission


TextbookOutcomes$READMISSION1[TextbookOutcomes$READMISSION1 == "NULL"] <- NA


TextbookOutcomes$READMISSION1[TextbookOutcomes$READMISSION1 == "NULL"] <- NA

table(TextbookOutcomes$READMISSION1)
## 
##    No   Yes 
## 27725  6219
## recoding -99 to NA's for length of stay. LOS considered >75th% other papers state median LOS for pancreatic cancer 7 days

TextbookOutcomes$DOptoDis[TextbookOutcomes$DOptoDis == -99] <- NA

table(TextbookOutcomes$PRNCPTX)
## 
## Distal Pancreatectomy             Lobectomy     Partial Lobectomy 
##                  8292                  4274                 11817 
##  Total Pancreatectomy      Trisegmentectomy               Whipple 
##                   706                  1366                 16464
TextbookOutcomes$RACE_NEW[TextbookOutcomes$RACE_NEW == "American Indian or Alaska Native"] <- "Other"
TextbookOutcomes$RACE_NEW[TextbookOutcomes$RACE_NEW == "Asian"] <- "Other"
TextbookOutcomes$RACE_NEW[TextbookOutcomes$RACE_NEW == "Native Hawaiian or Pacific Islander"] <- "Other"
TextbookOutcomes$RACE_NEW[TextbookOutcomes$RACE_NEW == "Black or African American"] <- "Black"
TextbookOutcomes$RACE_NEW[TextbookOutcomes$RACE_NEW == "Unknown/Not Reported"] <- NA
TextbookOutcomes$WEIGHT[TextbookOutcomes$WEIGHT ==-99] <- NA
TextbookOutcomes$HEIGHT[TextbookOutcomes$HEIGHT ==-99] <- NA
TextbookOutcomes$BMI <- (TextbookOutcomes$WEIGHT / TextbookOutcomes$HEIGHT / TextbookOutcomes$HEIGHT) *703
TextbookOutcomes <- 
  TextbookOutcomes %>%
  mutate(TextbookOutcomes, BMILabled = ifelse(BMI >= 40, "Morbidly Obese", 
                                               ifelse(BMI >= 30, "Obese", 
                                                      ifelse(BMI >= 25, "Overweight", 
                                                             ifelse(BMI >=18.5, "Normal Weight",
                                                                    ifelse(BMI <= 18.49, "Underweight", NA))))))
TextbookOutcomes$ASACLAS[TextbookOutcomes$ASACLAS == "5-Moribund"] <- NA

TextbookOutcomes$ASACLAS[TextbookOutcomes$ASACLAS == "None assigned"] <- NA
### Create MFI column

TextbookOutcomes$DIABETES[TextbookOutcomes$DIABETES == "INSULIN"] <- "Yes"

TextbookOutcomes$DIABETES[TextbookOutcomes$DIABETES == "NON-INSULIN"] <- "Yes"
###Create new numeric column that will be used to create the MFI5 column 0, >1, and greater than >2
##Diabeties recoded for MFI
TextbookOutcomes <-
  TextbookOutcomes %>%
  mutate(TextbookOutcomes, Diabetes1 = ifelse(DIABETES == "INSULIN", "Yes",
                                          ifelse(DIABETES == "NON-INSULIN", "Yes",
                                                 ifelse(DIABETES == "NO", "No", "NA"))))

###Create new numeric column that will be used to create the MFI5 column 0, >1, and greater than >2 ###CHF recoded for MFI

TextbookOutcomes <- TextbookOutcomes %>% mutate(TextbookOutcomes, CHFMFI = ifelse(HXCHF == “Yes”, “1”, ifelse(HXCHF == “No”, “0”, “NA”)))

###Create new numeric column that will be used to create the MFI5 column 0, >1, and greater than >2 ### COPD recoded for MFI

TextbookOutcomes <- TextbookOutcomes %>% mutate(TextbookOutcomes, COPDMFI = ifelse(HXCOPD == “Yes”, “1”, ifelse(HXCOPD == “No”, “0”, “NA”)))

###Create new numeric column that will be used to create the MFI5 column 0, >1, and greater than >2 ### Hypertension with medications recoded for MFI TextbookOutcomes <- TextbookOutcomes %>% mutate(TextbookOutcomes, HYPERMEDMFI = ifelse(HYPERMED == “Yes”, “1”, ifelse(HYPERMED == “No”, “0”, “NA”)))

##recode fnstatus2

TextbookOutcomes$FNSTATUS2[TextbookOutcomes$FNSTATUS2 == "Unknown"] <- NA

TextbookOutcomes <- TextbookOutcomes %>% mutate(TextbookOutcomes, FuncStatus = ifelse(FNSTATUS2 == “Yes”, “1”, ifelse(FNSTATUS2 == “No”, “0”, NA)))

##making the columns numeric to add them properly TextbookOutcomes\(HYPERMEDMFI <- as.numeric(TextbookOutcomes\)HYPERMEDMFI) TextbookOutcomes\(COPDMFI <- as.numeric(TextbookOutcomes\)COPDMFI) TextbookOutcomes\(CHFMFI <- as.numeric(TextbookOutcomes\)CHFMFI) TextbookOutcomes\(DiabetesMFI <- as.numeric(TextbookOutcomes\)DiabetesMFI) TextbookOutcomes\(FuncStatus <- as.numeric(TextbookOutcomes\)FuncStatus)

###Creating MFI column

TextbookOutcomes\(MFI5 <- (TextbookOutcomes\)HYPERMEDMFI + TextbookOutcomes\(COPDMFI + TextbookOutcomes\)CHFMFI + TextbookOutcomes\(DiabetesMFI + TextbookOutcomes\)FuncStatus )

##recoding Hep and Pan targeted Complications to yes or no (will need to add in targeted complications Pan: Fistula, Delayed Gastric Emptying and what about percdrain? | Hep: Bile leakage, hepliverfail)

TextbookOutcomes$PAN_FISTULA[TextbookOutcomes$PAN_FISTULA =="Biochemical Leak only"] <- "No"
TextbookOutcomes$PAN_FISTULA[TextbookOutcomes$PAN_FISTULA =="No evidence of Biochemical Leak or POPF"] <- "No"

TextbookOutcomes$PAN_FISTULA[TextbookOutcomes$PAN_FISTULA =="Yes-persistent drainage, drain continued >7 days"] <- "Yes"
TextbookOutcomes$PAN_FISTULA[TextbookOutcomes$PAN_FISTULA =="Yes-clinical diagnosis, percutaneous drainage performed"] <- "Yes"
TextbookOutcomes$PAN_FISTULA[TextbookOutcomes$PAN_FISTULA =="Yes-persistent drainage, percutaneous drainage performed"] <- "Yes"
TextbookOutcomes$PAN_FISTULA[TextbookOutcomes$PAN_FISTULA =="Yes-persistent drainage, NPO-TPN"] <- "Yes"
TextbookOutcomes$PAN_FISTULA[TextbookOutcomes$PAN_FISTULA =="Yes-clinical diagnosis, drain continued >7 days"] <- "Yes"
TextbookOutcomes$PAN_FISTULA[TextbookOutcomes$PAN_FISTULA =="Yes, Grade B POPF present"] <- "Yes"
TextbookOutcomes$PAN_FISTULA[TextbookOutcomes$PAN_FISTULA =="Yes-clinical diagnosis, spontaneous wound drainage"] <- "Yes"
TextbookOutcomes$PAN_FISTULA[TextbookOutcomes$PAN_FISTULA =="Yes-clinical diagnosis, NPO-TPN"] <- "Yes"

TextbookOutcomes$PAN_FISTULA[TextbookOutcomes$PAN_FISTULA =="Yes-clinical diagnosis, reoperation performed"] <- "Yes"
TextbookOutcomes$PAN_FISTULA[TextbookOutcomes$PAN_FISTULA =="Yes, Grade C POPF present"] <- "Yes"
TextbookOutcomes$PAN_FISTULA[TextbookOutcomes$PAN_FISTULA =="Yes-persistent drainage, reoperation performed"] <- "Yes"
TextbookOutcomes$PAN_FISTULA[TextbookOutcomes$PAN_FISTULA =="Unknown"] <- NA



####Delayed Gastric Emptying Recode to yes or no
TextbookOutcomes$PAN_DELGASTRIC_20140315[TextbookOutcomes$PAN_DELGASTRIC_20140315 =="Yes-no oral intake by POD 14"] <- "Yes"
TextbookOutcomes$PAN_DELGASTRIC_20140315[TextbookOutcomes$PAN_DELGASTRIC_20140315 =="Yes-tube to external drainage/NG tube present/reinserted"] <- "Yes"

TextbookOutcomes$PAN_DELGASTRIC_20140315[TextbookOutcomes$PAN_DELGASTRIC_20140315 =="Unknown"] <- NA

###recoding Hep targeted complications bileleakage
TextbookOutcomes$HEP_BILELEAKAGE[TextbookOutcomes$HEP_BILELEAKAGE =="Unknown"] <- NA
TextbookOutcomes$HEP_BILELEAKAGE[TextbookOutcomes$HEP_BILELEAKAGE =="Yes-clinical diagnosis, drain continued on or after POD3"] <- "Yes"
TextbookOutcomes$HEP_BILELEAKAGE[TextbookOutcomes$HEP_BILELEAKAGE =="Yes-clinical diagnosis, percutaneous drainage performed"] <- "Yes"
TextbookOutcomes$HEP_BILELEAKAGE[TextbookOutcomes$HEP_BILELEAKAGE =="Yes-clinical diagnosis, reoperation performed"] <- "Yes"
TextbookOutcomes$HEP_BILELEAKAGE[TextbookOutcomes$HEP_BILELEAKAGE =="Yes-clinical diagnosis, spontaneous wound drainage"] <- "Yes"
TextbookOutcomes$HEP_BILELEAKAGE[TextbookOutcomes$HEP_BILELEAKAGE =="Yes-persistent drainage, drain continued on or after POD3"] <- "Yes"
TextbookOutcomes$HEP_BILELEAKAGE[TextbookOutcomes$HEP_BILELEAKAGE =="Yes-persistent drainage, percutaneous drainage performed"] <- "Yes"
TextbookOutcomes$HEP_BILELEAKAGE[TextbookOutcomes$HEP_BILELEAKAGE =="Yes-persistent drainage, reoperation performed"] <- "Yes"

##recoding hep targeted complications Liver failure (Post Hepatectomy Liver Failure)
TextbookOutcomes$HEP_LIVERFAIL[TextbookOutcomes$HEP_LIVERFAIL =="Yes-PHLF (receiving clotting factors to maintain INR)"] <- "Yes"
TextbookOutcomes$HEP_LIVERFAIL[TextbookOutcomes$HEP_LIVERFAIL =="Yes-meets criteria for PHLF"] <- "Yes"

TextbookOutcomes$HEP_LIVERFAIL[TextbookOutcomes$HEP_LIVERFAIL =="No-does not meet criteria for PHLF"] <- "No"


table(TextbookOutcomes$HEP_LIVERFAIL)
## 
##    No   Yes 
## 16629   828
table(TextbookOutcomes$PRNCPTX)
## 
## Distal Pancreatectomy             Lobectomy     Partial Lobectomy 
##                  8292                  4274                 11817 
##  Total Pancreatectomy      Trisegmentectomy               Whipple 
##                   706                  1366                 16464
###Creating Targeted Column for Hep Complications
TextbookOutcomes1 <- TextbookOutcomes %>%
  mutate(TextbookOutcomes, HepComplications = ifelse(HEP_LIVERFAIL == "Yes", "Yes",
                                                   ifelse(HEP_BILELEAKAGE == "Yes", "Yes", "No")))
table(TextbookOutcomes1$HepComplications)
## 
##    No   Yes 
## 15510  1842
###Creating Targeted Column for Hep Complications
TextbookOutcomes2 <- TextbookOutcomes1 %>%
  mutate(TextbookOutcomes, PanComplications = ifelse(PAN_FISTULA == "Yes", "Yes",
                                                   ifelse(PAN_DELGASTRIC_20140315 == "Yes", "Yes", "No")))

str(TextbookOutcomes2$PanComplications)
##  chr [1:44235] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA ...
sum(is.na(TextbookOutcomes2$PanComplications))
## [1] 17754
###creating a true textbookoutcomes column for pan and liver

TextbookOutcomes2$PanComplications[TextbookOutcomes2$PanComplications ==NA] <- "NA"
TextbookOutcomes2$PanComplications[TextbookOutcomes2$HepComplications ==NA] <- "NA"


table(TextbookOutcomes2$PanComplications)
## 
##    No   Yes 
## 19879  6602
TextbookOutcomes2$PanComplications[TextbookOutcomes2$PanComplications ==NA] <- "NA"

sum(is.na(TextbookOutcomes2$PanComplications))
## [1] 17754
table(is.na(TextbookOutcomes2$PanComplications))
## 
## FALSE  TRUE 
## 26481 17754
### Ifelse function to create new column for any complication (will need to add in targeted complications Pan: Fistula, Delayed Gastric Emptying | Hep: Bile leakage, hepliverfail)
TextbookOutcomes3 <- TextbookOutcomes2 %>%
  mutate(TextbookOutcomes, AnyComplication = ifelse(CDARREST == "Cardiac Arrest Requiring CPR", "Yes",
                                                   ifelse(CDMI == "Myocardial Infarction", "Yes",
                                                          ifelse(OUPNEUMO == "Pneumonia", "Yes",
                                                                 ifelse(RENAINSF == "Progressive Renal Insufficiency", "Yes",
                                                                        ifelse(OPRENAFL == "Acute Renal Failure", "Yes",
                                                                               ifelse(PULEMBOL == "Pulmonary Embolism", "Yes",
                                                                                      ifelse(OTHDVT == "DVT Requiring Therapy", "Yes",
                                                                                             ifelse(RETURNOR == "Yes", "Yes",
                                                                                                    ifelse(ORGSPCSSI == "Organ/Space SSI", "Yes",
                                                                                                           ifelse(PRSEPIS == "Yes", "Yes",
                                                                                                                  ifelse(REINTUB == "Unplanned Intubation", "Yes",
                                                                                                                         ifelse(URNINFEC == "Urinary Tract Infection", "Yes",
                                                                                                                                ifelse(SUPINFEC == "Superficial Incisional SSI", "Yes",
                                                                                                                                       ifelse(FAILWEAN == "On Ventilator greater than 48 Hours", "Yes",
                                                                                                                                              ifelse(CNSCVA == "Stroke/CVA", "Yes",
                                                                                                                                                     ifelse(DEHIS == "Wound Disruption", "Yes", "No")))))))))))))))))
##subsetting Procedure with DoptoDIS information to obtain the mean

TextbookOutcomesDistalPan <-
  subset(TextbookOutcomes,  PRNCPTX == "Distal Pancreatectomy")

TextbookOutcomesLobectomy <-
  subset(TextbookOutcomes,  PRNCPTX == "Lobectomy")
  
TextbookOutcomesOtherPan <-
  subset(TextbookOutcomes,  PRNCPTX == "Other Pancreas")
  
TextbookOutcomesPartialLobectomy <-
  subset(TextbookOutcomes,  PRNCPTX == "Partial Lobectomy")
  
TextbookOutcomesTotalPan <-
  subset(TextbookOutcomes,  PRNCPTX == "Total Pancreatectomy")
  
TextbookOutcomesTrisegmentectomy <-
  subset(TextbookOutcomes,  PRNCPTX == "Trisegmentectomy")
  
TextbookOutcomesWhipple <-
  subset(TextbookOutcomes,  PRNCPTX == "Whipple")
###obtaining the 75th percential for LOS defined as day of surgery to date of discharge

SeventyFifthPercentileDistalPan = TextbookOutcomesDistalPan$DOptoDis
SeventyFifthPercentileLobectomy = TextbookOutcomesLobectomy$DOptoDis
SeventyFifthPercentileOtherPan = TextbookOutcomesOtherPan$DOptoDis
SeventyFifthPercentilePartialLobectomy = TextbookOutcomesPartialLobectomy$DOptoDis
SeventyFifthPercentileTotalPan = TextbookOutcomesTotalPan$DOptoDis
SeventyFifthPercentileTrisegmentectomy = TextbookOutcomesTrisegmentectomy$DOptoDis
SeventyFifthPercentileWhipple = TextbookOutcomesWhipple$DOptoDis


quantile(SeventyFifthPercentileDistalPan, na.rm = T)
##   0%  25%  50%  75% 100% 
##    0    4    5    7   88
quantile(SeventyFifthPercentileLobectomy, na.rm = T)
##   0%  25%  50%  75% 100% 
##    0    5    6    8  114
quantile(SeventyFifthPercentileOtherPan, na.rm = T)
##   0%  25%  50%  75% 100% 
##   NA   NA   NA   NA   NA
quantile(SeventyFifthPercentilePartialLobectomy, na.rm = T)
##   0%  25%  50%  75% 100% 
##    0    3    5    6   85
quantile(SeventyFifthPercentileTotalPan, na.rm = T)
##   0%  25%  50%  75% 100% 
##    0    6    8   12   89
quantile(SeventyFifthPercentileTrisegmentectomy, na.rm = T)
##   0%  25%  50%  75% 100% 
##    0    5    7   10   78
quantile(SeventyFifthPercentileWhipple, na.rm = T)
##   0%  25%  50%  75% 100% 
##    0    6    8   12  111
TextbookOutcomes3$ELOSDistalPan <- ifelse(TextbookOutcomes3$PRNCPTX == "Distal Pancreatectomy" & TextbookOutcomes3$DOptoDis >= 7, "Yes", "No")

TextbookOutcomes3$ELOSLobectomy <- ifelse(TextbookOutcomes3$PRNCPTX == "Lobectomy" & TextbookOutcomes3$DOptoDis >= 8, "Yes", "No")

TextbookOutcomes3$ELOSOtherPan <- ifelse(TextbookOutcomes3$PRNCPTX == "Other Pancreas" & TextbookOutcomes3$DOptoDis >= 7, "Yes", "No")

TextbookOutcomes3$ELOSPartial <- ifelse(TextbookOutcomes3$PRNCPTX == "Partial Lobectomy" & TextbookOutcomes3$DOptoDis >= 6, "Yes", "No")

TextbookOutcomes3$ELOSTotal <- ifelse(TextbookOutcomes3$PRNCPTX == "Total Pancreatectomy" & TextbookOutcomes3$DOptoDis >= 12, "Yes", "No")

TextbookOutcomes3$ELOSTrisegmentectomy <- ifelse(TextbookOutcomes3$PRNCPTX == "Trisegmentectomy" & TextbookOutcomes3$DOptoDis >= 10, "Yes", "No")

TextbookOutcomes3$ELOSWhipple <- ifelse(TextbookOutcomes3$PRNCPTX == "Whipple" & TextbookOutcomes3$DOptoDis >= 12, "Yes", "No")

table(TextbookOutcomes3$ELOSPartial)
## 
##    No   Yes 
## 39501  4153
###recode 9 and up to be considered increased LOS

PanTextbookOutcomes <- TextbookOutcomes3 %>%
  mutate(TextbookOutcomes3, ELOS = ifelse(ELOSDistalPan == "Yes", "Yes", 
                                          ifelse(ELOSLobectomy == "Yes", "Yes",
                                                 ifelse(ELOSOtherPan == "Yes", "Yes",
                                                        ifelse(ELOSPartial == "Yes", "Yes",
                                                               ifelse(ELOSTotal == "Yes", "Yes",
                                                                      ifelse(ELOSTrisegmentectomy == "Yes", "Yes",
                                                                             ifelse(ELOSWhipple == "Yes", "Yes", "No"))))))))
table(PanTextbookOutcomes$ELOS)
## 
##    No   Yes 
## 30578 12853

###subset pancreatic cancer patients prior to creating a Textbook Outcomes Table double the 75th percentile number afterwards

PanTextbookOutcomes <- subset(TextbookOutcomes4, PanComplications == “Yes” | PanComplications == “No”)

table(PanTextbookOutcomes$HepComplications)

PanTextbookOutcomes$TextbookOutcomeHep <- ifelse(PanTextbookOutcomes$ELOS == "No" & PanTextbookOutcomes$DOpertoD == "No" & PanTextbookOutcomes$HepComplications == "No" & PanTextbookOutcomes$AnyComplication == "No" & PanTextbookOutcomes$READMISSION1 == "No",
                  "Yes",
                  "No")



PanTextbookOutcomes$TextbookOutcomePan <- ifelse(PanTextbookOutcomes$ELOS == "No" & PanTextbookOutcomes$DOpertoD == "No" & PanTextbookOutcomes$PanComplications == "No" & PanTextbookOutcomes$AnyComplication == "No" & PanTextbookOutcomes$READMISSION1 == "No",
                  "Yes",
                  "No")

PanTextbookOutcomes$TextbookOutcomeAll <- ifelse(PanTextbookOutcomes$TextbookOutcomePan == "No" & PanTextbookOutcomes$TextbookOutcomeHep == "No",
                  "No",
                  "Yes")


table(PanTextbookOutcomes$TextbookOutcomeHep)
## 
##    No   Yes 
## 20105  6950
table(PanTextbookOutcomes$TextbookOutcomePan)
## 
##    No   Yes 
## 21045  9165
table(PanTextbookOutcomes$TextbookOutcomeAll)
## 
##    No   Yes 
## 19790 16115

PanTextbookOutcomes3 <- PanTextbookOutcomes %>% mutate(PanTextbookOutcomes, TextbookOutcomeAll = ifelse(TextbookOutcomePan == NA, “Yes”, ifelse(TextbookOutcomeHep == “Yes”, “Yes”,“NA”)))

table(PanTextbookOutcomes3$TextbookOutcomeAll)

PanTextbookOutcomes1 <- PanTextbookOutcomes %>% mutate(PanTextbookOutcomes, ModifiedFrailtyIndex5 = ifelse(MFI5 == 0, “0”, ifelse(MFI5 == 1, “1”, ifelse(MFI5 > 2, “>2”, NA))))

PantextBookOutcomesDT <-
  PanTextbookOutcomes %>%
  select(PRNCPTX, PatientGroup, Age, SEX, RACE_NEW, BMILabled, ASACLAS, DOpertoD, OPTIME, ELOS, READMISSION1, STEROID, DIABETES, HYPERMED, ASACLAS, ASCITES, HXCHF, DYSPNEA, SMOKE, HXCOPD, DIALYSIS, OPRENAFL, DISCANCR, FNSTATUS2, AnyComplication, TextbookOutcomeHep, TextbookOutcomePan, TextbookOutcomeAll, PanComplications, HepComplications)

table(PantextBookOutcomesDT$PatientGroup)
## 
##    Liver Pancreas 
##    17457    25462
PantextBookOutcomesDT %>%
  tbl_summary(
    by = PatientGroup,
    digits = all_continuous() ~ 2,) %>%
  add_p(pvalue_fun = ~style_pvalue(.x, digits = 2))%>%
  bold_p() %>%
  add_overall() %>%
  modify_header(label ~ "**Variable**") %>%
  modify_spanning_header(c("stat_1", "stat_2") ~ "**Textbook Outcomes*") %>%
  bold_labels() 
## 1316 observations missing `PatientGroup` have been removed. To include these observations, use `forcats::fct_explicit_na()` on `PatientGroup` column before passing to `tbl_summary()`.
## There was an error in 'add_p()' for variable 'PanComplications', p-value omitted:
## Error in stats::chisq.test(data[[variable]], as.factor(data[[by]])): 'x' and 'y' must have at least 2 levels
## There was an error in 'add_p()' for variable 'HepComplications', p-value omitted:
## Error in stats::chisq.test(data[[variable]], as.factor(data[[by]])): 'x' and 'y' must have at least 2 levels
Variable Overall, N = 42,9191 *Textbook Outcomes p-value2
Liver, N = 17,4571 Pancreas, N = 25,4621
PRNCPTX <0.001
Distal Pancreatectomy 8,292 (19%) 0 (0%) 8,292 (33%)
Lobectomy 4,274 (10.0%) 4,274 (24%) 0 (0%)
Partial Lobectomy 11,817 (28%) 11,817 (68%) 0 (0%)
Total Pancreatectomy 706 (1.6%) 0 (0%) 706 (2.8%)
Trisegmentectomy 1,366 (3.2%) 1,366 (7.8%) 0 (0%)
Whipple 16,464 (38%) 0 (0%) 16,464 (65%)
Age <0.001
0-65 22,630 (53%) 10,737 (62%) 11,893 (47%)
65-74 13,302 (31%) 4,660 (27%) 8,642 (34%)
75-84 6,399 (15%) 1,907 (11%) 4,492 (18%)
85+ 588 (1.4%) 153 (0.9%) 435 (1.7%)
SEX 0.040
female 21,264 (50%) 8,754 (50%) 12,510 (49%)
male 21,655 (50%) 8,703 (50%) 12,952 (51%)
RACE_NEW <0.001
Black 3,648 (10%) 1,468 (11%) 2,180 (9.8%)
Other 2,498 (6.9%) 1,298 (9.3%) 1,200 (5.4%)
White 29,943 (83%) 11,130 (80%) 18,813 (85%)
Unknown 6,830 3,561 3,269
BMILabled <0.001
Morbidly Obese 1,926 (4.5%) 908 (5.2%) 1,018 (4.0%)
Normal Weight 13,726 (32%) 5,229 (30%) 8,497 (34%)
Obese 11,570 (27%) 5,058 (29%) 6,512 (26%)
Overweight 14,594 (34%) 5,880 (34%) 8,714 (34%)
Underweight 879 (2.1%) 268 (1.5%) 611 (2.4%)
Unknown 224 114 110
ASACLAS <0.001
1-No Disturb 428 (1.0%) 253 (1.5%) 175 (0.7%)
2-Mild Disturb 10,403 (24%) 4,395 (25%) 6,008 (24%)
3-Severe Disturb 29,056 (68%) 11,495 (66%) 17,561 (69%)
4-Life Threat 2,951 (6.9%) 1,270 (7.3%) 1,681 (6.6%)
Unknown 81 44 37
DOpertoD 609 (1.4%) 257 (1.5%) 352 (1.4%) 0.47
OPTIME 270.00 (186.00, 373.00) 216.00 (154.00, 300.00) 312.00 (221.00, 409.00) <0.001
ELOS 12,853 (30%) 5,727 (33%) 7,126 (28%) <0.001
Unknown 242 50 192
READMISSION1 6,018 (18%) 1,725 (13%) 4,293 (22%) <0.001
Unknown 9,960 4,301 5,659
STEROID 1,290 (3.0%) 533 (3.1%) 757 (3.0%) 0.65
DIABETES <0.001
NO 33,197 (77%) 14,381 (82%) 18,816 (74%)
Yes 9,722 (23%) 3,076 (18%) 6,646 (26%)
HYPERMED 21,343 (50%) 8,021 (46%) 13,322 (52%) <0.001
ASCITES 169 (0.4%) 96 (0.5%) 73 (0.3%) <0.001
HXCHF 164 (0.4%) 59 (0.3%) 105 (0.4%) 0.25
DYSPNEA 0.082
AT REST 69 (0.2%) 37 (0.2%) 32 (0.1%)
MODERATE EXERTION 2,241 (5.2%) 921 (5.3%) 1,320 (5.2%)
No 40,609 (95%) 16,499 (95%) 24,110 (95%)
SMOKE 7,065 (16%) 2,698 (15%) 4,367 (17%) <0.001
HXCOPD 1,670 (3.9%) 640 (3.7%) 1,030 (4.0%) 0.049
DIALYSIS 151 (0.4%) 51 (0.3%) 100 (0.4%) 0.10
OPRENAFL 0.30
Acute Renal Failure 373 (0.9%) 162 (0.9%) 211 (0.8%)
No Complication 42,546 (99%) 17,295 (99%) 25,251 (99%)
DISCANCR 8,465 (20%) 7,047 (40%) 1,418 (5.6%) <0.001
FNSTATUS2 0.077
Independent 42,512 (99%) 17,308 (99%) 25,204 (99%)
Partially Dependent 319 (0.7%) 111 (0.6%) 208 (0.8%)
Totally Dependent 24 (<0.1%) 8 (<0.1%) 16 (<0.1%)
Unknown 64 30 34
AnyComplication 10,522 (25%) 3,269 (19%) 7,253 (28%) <0.001
TextbookOutcomeHep 6,950 (26%) 6,950 (47%) 0 (0%) <0.001
Unknown 16,209 2,801 13,408
TextbookOutcomePan 8,759 (30%) 0 (0%) 8,759 (40%) <0.001
Unknown 13,555 10,066 3,489
TextbookOutcomeAll 15,709 (45%) 6,950 (48%) 8,759 (42%) <0.001
Unknown 7,765 3,116 4,649
PanComplications 6,370 (25%) 0 (NA%) 6,370 (25%)
Unknown 17,737 17,457 280
HepComplications 1,842 (11%) 1,842 (11%) 0 (NA%)
Unknown 25,567 105 25,462

1 Statistics presented: n (%); Median (IQR)

2 Statistical tests performed: chi-square test of independence; Wilcoxon rank-sum test 

TextbookOutcomeOR <- 
  PantextBookOutcomesDT %>%
  select(PRNCPTX, PatientGroup, Age, SEX, RACE_NEW, BMILabled, ASACLAS, DOpertoD, OPTIME, ELOS, READMISSION1, STEROID, DIABETES, HYPERMED, ASACLAS, ASCITES, HXCHF, DYSPNEA, SMOKE, HXCOPD, DIALYSIS, OPRENAFL, DISCANCR, FNSTATUS2, AnyComplication, TextbookOutcomeHep, TextbookOutcomePan, TextbookOutcomeAll, PanComplications, HepComplications)
TextbookOutcomeOR <- TextbookOutcomeOR %>%
  mutate(TextbookOutcomeOR, TextbookOutcomePanOR = ifelse(TextbookOutcomePan =="Yes", "1",
                                                  ifelse(TextbookOutcomePan =="No", "0", NA)))

TextbookOutcomeOR$TextbookOutcomePanOR <- as.numeric(TextbookOutcomeOR$TextbookOutcomePanOR)


TextbookOutcomeOR <- TextbookOutcomeOR %>%
  mutate(TextbookOutcomeOR, TextbookOutcomeHepOR = ifelse(TextbookOutcomeHep =="Yes", "1",
                                                  ifelse(TextbookOutcomeHep =="No", "0", NA)))

TextbookOutcomeOR <- TextbookOutcomeOR %>%
  mutate(TextbookOutcomeOR, TextbookOutcomeAllOR = ifelse(TextbookOutcomeAll =="Yes", "1",
                                                  ifelse(TextbookOutcomeAll =="No", "0", NA)))

TextbookOutcomeOR$TextbookOutcomeHepOR <- as.numeric(TextbookOutcomeOR$TextbookOutcomeHepOR)
TextbookOutcomeOR$TextbookOutcomeAllOR <- as.numeric(TextbookOutcomeOR$TextbookOutcomeAllOR)
TextbookOutcomeOR$TextbookOutcomePanOR <- as.numeric(TextbookOutcomeOR$TextbookOutcomePanOR)
PanTextbookOutcomesSubset <- 
  subset(TextbookOutcomeOR, PatientGroup == "Pancreas")
HepTextbookOutcomesSubset <- 
  subset(TextbookOutcomeOR, PatientGroup == "Liver")

table(PanTextbookOutcomesSubset$TextbookOutcomePanOR)
## 
##     0     1 
## 13214  8759
###changling the referent level

PanTextbookOutcomesSubset$PRNCPTX <- as.factor(PanTextbookOutcomesSubset$PRNCPTX)
PanTextbookOutcomesSubset$SEX <- as.factor(PanTextbookOutcomesSubset$SEX)
PanTextbookOutcomesSubset$ASACLAS <- as.factor(PanTextbookOutcomesSubset$ASACLAS)
PanTextbookOutcomesSubset$BMILabled <- as.factor(PanTextbookOutcomesSubset$BMILabled)
PanTextbookOutcomesSubset$STEROID <- as.factor(PanTextbookOutcomesSubset$STEROID)
PanTextbookOutcomesSubset$DIABETES <- as.factor(PanTextbookOutcomesSubset$DIABETES)
PanTextbookOutcomesSubset$HYPERMED <- as.factor(PanTextbookOutcomesSubset$HYPERMED)
PanTextbookOutcomesSubset$ASCITES <- as.factor(PanTextbookOutcomesSubset$ASCITES)
PanTextbookOutcomesSubset$HXCOPD <- as.factor(PanTextbookOutcomesSubset$HXCOPD)
PanTextbookOutcomesSubset$DIALYSIS <- as.factor(PanTextbookOutcomesSubset$DIALYSIS)
PanTextbookOutcomesSubset$DISCANCR <- as.factor(PanTextbookOutcomesSubset$DISCANCR)
PanTextbookOutcomesSubset$HXCHF <- as.factor(PanTextbookOutcomesSubset$HXCHF)
PanTextbookOutcomesSubset$Age <- as.factor(PanTextbookOutcomesSubset$Age)
PanTextbookOutcomesSubset$FNSTATUS2 <- as.factor(PanTextbookOutcomesSubset$FNSTATUS2)
PanTextbookOutcomesSubset$SMOKE <- as.factor(PanTextbookOutcomesSubset$SMOKE)
PanTextbookOutcomesSubset$PatientGroup <- as.factor(PanTextbookOutcomesSubset$PatientGroup)


PanTextbookOutcomesSubset$PRNCPTX = relevel(PanTextbookOutcomesSubset$PRNCPTX, ref="Distal Pancreatectomy")
PanTextbookOutcomesSubset$SEX = relevel(PanTextbookOutcomesSubset$SEX, ref="male")
PanTextbookOutcomesSubset$BMILabled = relevel(PanTextbookOutcomesSubset$BMILabled, ref="Overweight")
PanTextbookOutcomesSubset$ASACLAS = relevel(PanTextbookOutcomesSubset$ASACLAS, ref="4-Life Threat")
PanTextbookOutcomesSubset$DIABETES = relevel(PanTextbookOutcomesSubset$DIABETES, ref="NO")
PanTextbookOutcomesSubset$HYPERMED = relevel(PanTextbookOutcomesSubset$HYPERMED, ref="Yes")
PanTextbookOutcomesSubset$STEROID = relevel(PanTextbookOutcomesSubset$STEROID, ref="Yes")
PanTextbookOutcomesSubset$ASCITES = relevel(PanTextbookOutcomesSubset$ASCITES, ref="Yes")
PanTextbookOutcomesSubset$HXCOPD = relevel(PanTextbookOutcomesSubset$HXCOPD, ref="Yes")
PanTextbookOutcomesSubset$DIALYSIS = relevel(PanTextbookOutcomesSubset$DIALYSIS, ref="Yes")
PanTextbookOutcomesSubset$SMOKE = relevel(PanTextbookOutcomesSubset$SMOKE, ref="No")
PanTextbookOutcomesSubset$DISCANCR = relevel(PanTextbookOutcomesSubset$DISCANCR, ref="Yes")
PanTextbookOutcomesSubset$HXCHF = relevel(PanTextbookOutcomesSubset$HXCHF, ref="Yes")
PanTextbookOutcomesSubset$Age = relevel(PanTextbookOutcomesSubset$Age, ref="85+")
PanTextbookOutcomesSubset$FNSTATUS2 = relevel(PanTextbookOutcomesSubset$FNSTATUS2, ref="Totally Dependent")
UVPan <- PanTextbookOutcomesSubset %>%
  select(PRNCPTX, Age, SEX, RACE_NEW, BMILabled, ASACLAS, STEROID, ASCITES, DYSPNEA, SMOKE, DIALYSIS, DISCANCR, DIABETES, HYPERMED, HXCOPD, HXCHF, TextbookOutcomePanOR) %>%
  tbl_uvregression(
    method = glm,
    y = TextbookOutcomePanOR,
    method.args = list(family = binomial),
    exponentiate = TRUE,
    pvalue_fun = ~style_pvalue(.x, digits = 2)
  ) %>%
  bold_p(t = 0.10) %>%
  bold_labels() %>%
  italicize_levels()
###changling the referent level

PanTextbookOutcomesSubset$PRNCPTX <- as.factor(PanTextbookOutcomesSubset$PRNCPTX)
PanTextbookOutcomesSubset$SEX <- as.factor(PanTextbookOutcomesSubset$SEX)
PanTextbookOutcomesSubset$ASACLAS <- as.factor(PanTextbookOutcomesSubset$ASACLAS)
PanTextbookOutcomesSubset$BMILabled <- as.factor(PanTextbookOutcomesSubset$BMILabled)
PanTextbookOutcomesSubset$STEROID <- as.factor(PanTextbookOutcomesSubset$STEROID)
PanTextbookOutcomesSubset$DIABETES <- as.factor(PanTextbookOutcomesSubset$DIABETES)
PanTextbookOutcomesSubset$HYPERMED <- as.factor(PanTextbookOutcomesSubset$HYPERMED)
PanTextbookOutcomesSubset$ASCITES <- as.factor(PanTextbookOutcomesSubset$ASCITES)
PanTextbookOutcomesSubset$HXCOPD <- as.factor(PanTextbookOutcomesSubset$HXCOPD)
PanTextbookOutcomesSubset$DIALYSIS <- as.factor(PanTextbookOutcomesSubset$DIALYSIS)
PanTextbookOutcomesSubset$DISCANCR <- as.factor(PanTextbookOutcomesSubset$DISCANCR)
PanTextbookOutcomesSubset$HXCHF <- as.factor(PanTextbookOutcomesSubset$HXCHF)
PanTextbookOutcomesSubset$Age <- as.factor(PanTextbookOutcomesSubset$Age)
PanTextbookOutcomesSubset$FNSTATUS2 <- as.factor(PanTextbookOutcomesSubset$FNSTATUS2)
PanTextbookOutcomesSubset$PatientGroup <- as.factor(PanTextbookOutcomesSubset$PatientGroup)


PanTextbookOutcomesSubset$PRNCPTX = relevel(PanTextbookOutcomesSubset$PRNCPTX, ref="Distal Pancreatectomy")
PanTextbookOutcomesSubset$SEX = relevel(PanTextbookOutcomesSubset$SEX, ref="male")
PanTextbookOutcomesSubset$ASACLAS = relevel(PanTextbookOutcomesSubset$ASACLAS, ref="4-Life Threat")
PanTextbookOutcomesSubset$BMILabled = relevel(PanTextbookOutcomesSubset$BMILabled, ref="Normal Weight")
PanTextbookOutcomesSubset$DIABETES = relevel(PanTextbookOutcomesSubset$DIABETES, ref="NO")
PanTextbookOutcomesSubset$HYPERMED = relevel(PanTextbookOutcomesSubset$HYPERMED, ref="Yes")
PanTextbookOutcomesSubset$STEROID = relevel(PanTextbookOutcomesSubset$STEROID, ref="Yes")
PanTextbookOutcomesSubset$ASCITES = relevel(PanTextbookOutcomesSubset$ASCITES, ref="Yes")
PanTextbookOutcomesSubset$HXCOPD = relevel(PanTextbookOutcomesSubset$HXCOPD, ref="Yes")
PanTextbookOutcomesSubset$DIALYSIS = relevel(PanTextbookOutcomesSubset$DIALYSIS, ref="Yes")
PanTextbookOutcomesSubset$DISCANCR = relevel(PanTextbookOutcomesSubset$DISCANCR, ref="Yes")
PanTextbookOutcomesSubset$HXCHF = relevel(PanTextbookOutcomesSubset$HXCHF, ref="Yes")
PanTextbookOutcomesSubset$Age = relevel(PanTextbookOutcomesSubset$Age, ref="85+")
PanTextbookOutcomesSubset$FNSTATUS2 = relevel(PanTextbookOutcomesSubset$FNSTATUS2, ref="Totally Dependent")
####Multivariant analysis serious complications
MVORPan <- 
  glm(TextbookOutcomePanOR ~ PRNCPTX + Age + SEX + RACE_NEW + BMILabled + ASACLAS + STEROID + ASCITES + DYSPNEA + SMOKE + DIALYSIS + DISCANCR + DIABETES + HYPERMED + HXCOPD + HXCHF, 
      data = PanTextbookOutcomesSubset,
      family = binomial("logit"), 
      na.action =na.omit)
###continue on to merge the table sets
MVORTablePan <-
  tbl_regression(MVORPan, exponentiate = T,
                 pvalue_fun = ~style_pvalue(.x, digits = 2),
  ) %>%
  bold_p(t = 0.10) %>%
  bold_labels() %>%
  italicize_levels() 
tbl_merge(
  list(UVPan, MVORTablePan),
  tab_spanner = c("**Univariable**", "**Multivariable**")
) %>%
  bold_labels() %>%
  italicize_levels()
Characteristic Univariable Multivariable
N OR1 95% CI1 p-value OR1 95% CI1 p-value
PRNCPTX 21,973
Distal Pancreatectomy — — — —
Total Pancreatectomy 1.07 0.90, 1.26 0.45 1.01 0.84, 1.21 0.90
Whipple 1.03 0.97, 1.09 0.36 1.04 0.97, 1.11 0.25
Age 21,973
85+ — — — —
0-65 1.48 1.19, 1.84 <0.001 1.49 1.19, 1.88 <0.001
65-74 1.40 1.13, 1.75 0.003 1.44 1.15, 1.82 0.002
75-84 1.22 0.98, 1.54 0.075 1.26 1.00, 1.60 0.052
SEX 21,973
male — — — —
female 1.21 1.15, 1.28 <0.001 1.20 1.13, 1.28 <0.001
RACE_NEW 19,104
Black — — — —
Other 0.83 0.70, 0.97 0.017 0.72 0.61, 0.84 <0.001
White 1.09 0.99, 1.20 0.082 1.09 0.98, 1.20 0.10
BMILabled 21,881
Overweight — — 0.87 0.81, 0.93 <0.001
Morbidly Obese 0.70 0.61, 0.81 <0.001 0.60 0.51, 0.70 <0.001
Normal Weight 1.16 1.09, 1.24 <0.001 — —
Obese 0.85 0.80, 0.92 <0.001 0.72 0.67, 0.78 <0.001
Underweight 1.27 1.06, 1.52 0.008 1.12 0.92, 1.35 0.25
ASACLAS 21,937
4-Life Threat — — — —
1-No Disturb 2.17 1.55, 3.03 <0.001 1.86 1.25, 2.78 0.002
2-Mild Disturb 1.93 1.71, 2.19 <0.001 1.75 1.49, 2.04 <0.001
3-Severe Disturb 1.53 1.36, 1.72 <0.001 1.47 1.27, 1.70 <0.001
STEROID 21,973
Yes — — — —
No 1.28 1.09, 1.52 0.003 1.22 1.02, 1.46 0.027
ASCITES 21,973
Yes — — — —
No 1.55 0.94, 2.64 0.094 1.56 0.91, 2.79 0.12
DYSPNEA 21,973
AT REST — — — —
MODERATE EXERTION 0.94 0.43, 2.20 0.87 1.00 0.42, 2.66 >0.99
No 1.35 0.62, 3.16 0.46 1.21 0.51, 3.20 0.68
SMOKE 21,973
No — — — —
Yes 1.04 0.97, 1.12 0.25 1.02 0.94, 1.10 0.68
DIALYSIS 21,973
Yes — — — —
No 2.55 1.55, 4.41 <0.001 2.09 1.21, 3.85 0.012
DISCANCR 21,973
Yes — — — —
No 1.54 1.37, 1.75 <0.001 1.49 1.30, 1.70 <0.001
DIABETES 21,973
NO — — — —
Yes 0.94 0.88, 1.00 0.047 1.08 1.00, 1.16 0.036
HYPERMED 21,973
Yes — — — —
No 1.19 1.13, 1.26 <0.001 1.04 0.98, 1.11 0.19
HXCOPD 21,973
Yes — — — —
No 1.62 1.40, 1.87 <0.001 1.47 1.26, 1.73 <0.001
HXCHF 21,973
Yes — — — —
No 2.38 1.49, 3.95 <0.001 1.56 0.94, 2.71 0.10

1 OR = Odds Ratio, CI = Confidence Interval 

###changling the referent level


HepTextbookOutcomesSubset$PRNCPTX <- as.factor(HepTextbookOutcomesSubset$PRNCPTX)
HepTextbookOutcomesSubset$SEX <- as.factor(HepTextbookOutcomesSubset$SEX)
HepTextbookOutcomesSubset$BMILabled <- as.factor(HepTextbookOutcomesSubset$BMILabled)
HepTextbookOutcomesSubset$ASACLAS <- as.factor(HepTextbookOutcomesSubset$ASACLAS)
HepTextbookOutcomesSubset$STEROID <- as.factor(HepTextbookOutcomesSubset$STEROID)
HepTextbookOutcomesSubset$DIABETES <- as.factor(HepTextbookOutcomesSubset$DIABETES)
HepTextbookOutcomesSubset$HYPERMED <- as.factor(HepTextbookOutcomesSubset$HYPERMED)
HepTextbookOutcomesSubset$ASCITES <- as.factor(HepTextbookOutcomesSubset$ASCITES)
HepTextbookOutcomesSubset$HXCOPD <- as.factor(HepTextbookOutcomesSubset$HXCOPD)
HepTextbookOutcomesSubset$DIALYSIS <- as.factor(HepTextbookOutcomesSubset$DIALYSIS)
HepTextbookOutcomesSubset$DISCANCR <- as.factor(HepTextbookOutcomesSubset$DISCANCR)
HepTextbookOutcomesSubset$HXCHF <- as.factor(HepTextbookOutcomesSubset$HXCHF)
HepTextbookOutcomesSubset$SMOKE <- as.factor(HepTextbookOutcomesSubset$SMOKE)
HepTextbookOutcomesSubset$Age <- as.factor(HepTextbookOutcomesSubset$Age)


HepTextbookOutcomesSubset$PRNCPTX = relevel(HepTextbookOutcomesSubset$PRNCPTX, ref="Trisegmentectomy")
HepTextbookOutcomesSubset$SEX = relevel(HepTextbookOutcomesSubset$SEX, ref="male")
HepTextbookOutcomesSubset$BMILabled = relevel(HepTextbookOutcomesSubset$BMILabled, ref="Normal Weight")
HepTextbookOutcomesSubset$ASACLAS = relevel(HepTextbookOutcomesSubset$ASACLAS, ref="4-Life Threat")
HepTextbookOutcomesSubset$DIABETES = relevel(HepTextbookOutcomesSubset$DIABETES, ref="Yes")
HepTextbookOutcomesSubset$HYPERMED = relevel(HepTextbookOutcomesSubset$HYPERMED, ref="Yes")
HepTextbookOutcomesSubset$STEROID = relevel(HepTextbookOutcomesSubset$STEROID, ref="Yes")
HepTextbookOutcomesSubset$ASCITES = relevel(HepTextbookOutcomesSubset$ASCITES, ref="Yes")
HepTextbookOutcomesSubset$HXCOPD = relevel(HepTextbookOutcomesSubset$HXCOPD, ref="Yes")
HepTextbookOutcomesSubset$DIALYSIS = relevel(HepTextbookOutcomesSubset$DIALYSIS, ref="Yes")
HepTextbookOutcomesSubset$DISCANCR = relevel(HepTextbookOutcomesSubset$DISCANCR, ref="Yes")
HepTextbookOutcomesSubset$HXCHF = relevel(HepTextbookOutcomesSubset$HXCHF, ref="Yes")
HepTextbookOutcomesSubset$SMOKE = relevel(HepTextbookOutcomesSubset$SMOKE, ref="Yes")
HepTextbookOutcomesSubset$Age = relevel(HepTextbookOutcomesSubset$Age, ref="85+")
UVHep <- HepTextbookOutcomesSubset %>%
  select(PRNCPTX, Age, SEX, RACE_NEW, BMILabled, ASACLAS, STEROID, ASCITES, DYSPNEA, SMOKE, DIALYSIS, DISCANCR, DIABETES, HYPERMED, HXCOPD, HXCHF, TextbookOutcomeHepOR) %>%
  tbl_uvregression(
    method = glm,
    y = TextbookOutcomeHepOR,
    method.args = list(family = binomial),
    exponentiate = TRUE,
    pvalue_fun = ~style_pvalue(.x, digits = 2)
  ) %>%
  bold_p(t = 0.10) %>%
  bold_labels() %>%
  italicize_levels()
HepTextbookOutcomesSubset$PRNCPTX <- as.factor(HepTextbookOutcomesSubset$PRNCPTX)
HepTextbookOutcomesSubset$BMILabled <- as.factor(HepTextbookOutcomesSubset$BMILabled)
HepTextbookOutcomesSubset$ASACLAS <- as.factor(HepTextbookOutcomesSubset$ASACLAS)
HepTextbookOutcomesSubset$STEROID <- as.factor(HepTextbookOutcomesSubset$STEROID)
HepTextbookOutcomesSubset$DIABETES <- as.factor(HepTextbookOutcomesSubset$DIABETES)
HepTextbookOutcomesSubset$HYPERMED <- as.factor(HepTextbookOutcomesSubset$HYPERMED)
HepTextbookOutcomesSubset$ASCITES <- as.factor(HepTextbookOutcomesSubset$ASCITES)
HepTextbookOutcomesSubset$HXCOPD <- as.factor(HepTextbookOutcomesSubset$HXCOPD)
HepTextbookOutcomesSubset$DIALYSIS <- as.factor(HepTextbookOutcomesSubset$DIALYSIS)
HepTextbookOutcomesSubset$DISCANCR <- as.factor(HepTextbookOutcomesSubset$DISCANCR)
HepTextbookOutcomesSubset$HXCHF <- as.factor(HepTextbookOutcomesSubset$HXCHF)
HepTextbookOutcomesSubset$SMOKE <- as.factor(HepTextbookOutcomesSubset$SMOKE)
HepTextbookOutcomesSubset$Age <- as.factor(HepTextbookOutcomesSubset$Age)


HepTextbookOutcomesSubset$PRNCPTX = relevel(HepTextbookOutcomesSubset$PRNCPTX, ref="Trisegmentectomy")
HepTextbookOutcomesSubset$BMILabled = relevel(HepTextbookOutcomesSubset$BMILabled, ref="Normal Weight")
HepTextbookOutcomesSubset$ASACLAS = relevel(HepTextbookOutcomesSubset$ASACLAS, ref="4-Life Threat")
HepTextbookOutcomesSubset$DIABETES = relevel(HepTextbookOutcomesSubset$DIABETES, ref="Yes")
HepTextbookOutcomesSubset$HYPERMED = relevel(HepTextbookOutcomesSubset$HYPERMED, ref="Yes")
HepTextbookOutcomesSubset$STEROID = relevel(HepTextbookOutcomesSubset$STEROID, ref="Yes")
HepTextbookOutcomesSubset$ASCITES = relevel(HepTextbookOutcomesSubset$ASCITES, ref="Yes")
HepTextbookOutcomesSubset$HXCOPD = relevel(HepTextbookOutcomesSubset$HXCOPD, ref="Yes")
HepTextbookOutcomesSubset$DIALYSIS = relevel(HepTextbookOutcomesSubset$DIALYSIS, ref="Yes")
HepTextbookOutcomesSubset$DISCANCR = relevel(HepTextbookOutcomesSubset$DISCANCR, ref="Yes")
HepTextbookOutcomesSubset$HXCHF = relevel(HepTextbookOutcomesSubset$HXCHF, ref="Yes")
HepTextbookOutcomesSubset$SMOKE = relevel(HepTextbookOutcomesSubset$SMOKE, ref="Yes")
HepTextbookOutcomesSubset$Age = relevel(HepTextbookOutcomesSubset$Age, ref="85+")
####Multivariant analysis serious complications
MVORHep <- 
  glm(TextbookOutcomeHepOR ~ PRNCPTX + Age + SEX + RACE_NEW + BMILabled + ASACLAS + STEROID + ASCITES + DYSPNEA + SMOKE + DIALYSIS + DISCANCR + DIABETES + HYPERMED + HXCOPD + HXCHF, 
      data = HepTextbookOutcomesSubset,
      family = binomial("logit"), 
      na.action =na.omit)
###continue on to merge the table sets
MVORTableHep <-
  tbl_regression(MVORHep, exponentiate = T,
                 pvalue_fun = ~style_pvalue(.x, digits = 2),
  ) %>%
  bold_p(t = 0.10) %>%
  bold_labels() %>%
  italicize_levels() 
tbl_merge(
  list(UVHep, MVORTableHep),
  tab_spanner = c("**Univariable**", "**Multivariable**")
) %>%
  bold_labels() %>%
  italicize_levels()
Characteristic Univariable Multivariable
N OR1 95% CI1 p-value OR1 95% CI1 p-value
PRNCPTX 14,656
Trisegmentectomy — — — —
Lobectomy 1.30 1.14, 1.49 <0.001 1.28 1.10, 1.50 0.002
Partial Lobectomy 1.36 1.20, 1.54 <0.001 1.37 1.19, 1.58 <0.001
Age 14,656
85+ — — — —
0-65 1.94 1.36, 2.82 <0.001 1.88 1.25, 2.88 0.003
65-74 1.40 0.98, 2.04 0.071 1.50 0.99, 2.30 0.057
75-84 1.09 0.75, 1.60 0.65 1.18 0.77, 1.82 0.46
SEX 14,656
male — — — —
female 1.36 1.27, 1.45 <0.001 1.26 1.16, 1.36 <0.001
RACE_NEW 11,659
Black — — — —
Other 0.96 0.82, 1.14 0.67 0.94 0.79, 1.12 0.47
White 1.11 0.99, 1.25 0.084 1.09 0.96, 1.23 0.17
BMILabled 14,563
Normal Weight — — — —
Morbidly Obese 0.96 0.82, 1.12 0.59 0.95 0.80, 1.13 0.55
Obese 1.03 0.95, 1.12 0.51 1.08 0.98, 1.20 0.12
Overweight 1.08 0.99, 1.17 0.068 1.15 1.04, 1.27 0.004
Underweight 0.79 0.60, 1.03 0.082 0.74 0.55, 1.00 0.049
ASACLAS 14,617
4-Life Threat — — — —
1-No Disturb 3.24 2.40, 4.39 <0.001 1.99 1.24, 3.25 0.005
2-Mild Disturb 2.62 2.28, 3.02 <0.001 1.81 1.49, 2.19 <0.001
3-Severe Disturb 1.76 1.55, 2.01 <0.001 1.38 1.15, 1.65 <0.001
STEROID 14,656
Yes — — — —
No 1.27 1.05, 1.54 0.013 1.16 0.94, 1.43 0.17
ASCITES 14,656
Yes — — — —
No 3.34 2.06, 5.71 <0.001 3.01 1.78, 5.37 <0.001
DYSPNEA 14,656
AT REST — — — —
MODERATE EXERTION 1.36 0.67, 2.93 0.40 1.22 0.58, 2.69 0.62
No 2.01 1.01, 4.28 0.055 1.49 0.72, 3.26 0.30
SMOKE 14,656
Yes — — — —
No 1.10 1.01, 1.20 0.038 1.05 0.95, 1.17 0.35
DIALYSIS 14,656
Yes — — — —
No 2.55 1.32, 5.31 0.008 1.78 0.87, 3.93 0.13
DISCANCR 14,656
Yes — — — —
No 1.14 1.07, 1.22 <0.001 1.17 1.08, 1.26 <0.001
DIABETES 14,656
Yes — — — —
NO 1.46 1.34, 1.59 <0.001 1.24 1.12, 1.38 <0.001
HYPERMED 14,656
Yes — — — —
No 1.34 1.26, 1.43 <0.001 1.10 1.01, 1.20 0.023
HXCOPD 14,656
Yes — — — —
No 1.64 1.38, 1.96 <0.001 1.25 1.02, 1.53 0.030
HXCHF 14,656
Yes — — — —
No 2.85 1.57, 5.55 <0.001 2.00 1.03, 4.17 0.050

1 OR = Odds Ratio, CI = Confidence Interval 

###changling the referent level

TextbookOutcomeOR$PatientGroup <- as.factor(TextbookOutcomeOR$PatientGroup)
TextbookOutcomeOR$BMILabled <- as.factor(TextbookOutcomeOR$BMILabled)
TextbookOutcomeOR$ASACLAS <- as.factor(TextbookOutcomeOR$ASACLAS)
TextbookOutcomeOR$STEROID <- as.factor(TextbookOutcomeOR$STEROID)
TextbookOutcomeOR$DIABETES <- as.factor(TextbookOutcomeOR$DIABETES)
TextbookOutcomeOR$HYPERMED <- as.factor(TextbookOutcomeOR$HYPERMED)
TextbookOutcomeOR$ASCITES <- as.factor(TextbookOutcomeOR$ASCITES)
TextbookOutcomeOR$HXCOPD <- as.factor(TextbookOutcomeOR$HXCOPD)
TextbookOutcomeOR$DIALYSIS <- as.factor(TextbookOutcomeOR$DIALYSIS)
TextbookOutcomeOR$DISCANCR <- as.factor(TextbookOutcomeOR$DISCANCR)
TextbookOutcomeOR$HXCHF <- as.factor(TextbookOutcomeOR$HXCHF)
TextbookOutcomeOR$SMOKE <- as.factor(TextbookOutcomeOR$SMOKE)
TextbookOutcomeOR$Age <- as.factor(TextbookOutcomeOR$Age)
TextbookOutcomeOR$SEX <- as.factor(TextbookOutcomeOR$SEX)




TextbookOutcomeOR$BMILabled = relevel(TextbookOutcomeOR$BMILabled, ref="Normal Weight")
TextbookOutcomeOR$ASACLAS = relevel(TextbookOutcomeOR$ASACLAS, ref="4-Life Threat")
TextbookOutcomeOR$DIABETES = relevel(TextbookOutcomeOR$DIABETES, ref="Yes")
TextbookOutcomeOR$SEX = relevel(TextbookOutcomeOR$SEX, ref="male")
TextbookOutcomeOR$SMOKE = relevel(TextbookOutcomeOR$SMOKE, ref="Yes")
TextbookOutcomeOR$STEROID = relevel(TextbookOutcomeOR$STEROID, ref="Yes")
TextbookOutcomeOR$ASCITES = relevel(TextbookOutcomeOR$ASCITES, ref="Yes")
TextbookOutcomeOR$HXCOPD = relevel(TextbookOutcomeOR$HXCOPD, ref="Yes")
TextbookOutcomeOR$DIALYSIS = relevel(TextbookOutcomeOR$DIALYSIS, ref="Yes")
TextbookOutcomeOR$DISCANCR = relevel(TextbookOutcomeOR$DISCANCR, ref="Yes")
TextbookOutcomeOR$HXCHF = relevel(TextbookOutcomeOR$HXCHF, ref="Yes")
TextbookOutcomeOR$Age = relevel(TextbookOutcomeOR$Age, ref="85+")
###Univariant, multivariant for everyone


UVAll <- TextbookOutcomeOR %>%
  select(PatientGroup, Age, SEX, RACE_NEW, BMILabled, ASACLAS, STEROID, ASCITES, DYSPNEA, SMOKE, DIALYSIS, DISCANCR, DIABETES, HYPERMED, HXCOPD, HXCHF, TextbookOutcomeAllOR) %>%
  tbl_uvregression(
    method = glm,
    y = TextbookOutcomeAllOR,
    method.args = list(family = binomial),
    exponentiate = TRUE,
    pvalue_fun = ~style_pvalue(.x, digits = 2)
  ) %>%
  bold_p(t = 0.10) %>%
  bold_labels() %>%
  italicize_levels()
###changling the referent level

TextbookOutcomeOR$PatientGroup <- as.factor(TextbookOutcomeOR$PatientGroup)
TextbookOutcomeOR$SEX <- as.factor(TextbookOutcomeOR$SEX)
TextbookOutcomeOR$ASACLAS <- as.factor(TextbookOutcomeOR$ASACLAS)
TextbookOutcomeOR$STEROID <- as.factor(TextbookOutcomeOR$STEROID)
TextbookOutcomeOR$DIABETES <- as.factor(TextbookOutcomeOR$DIABETES)
TextbookOutcomeOR$SMOKE <- as.factor(TextbookOutcomeOR$SMOKE)
TextbookOutcomeOR$ASCITES <- as.factor(TextbookOutcomeOR$ASCITES)
TextbookOutcomeOR$HXCOPD <- as.factor(TextbookOutcomeOR$HXCOPD)
TextbookOutcomeOR$DIALYSIS <- as.factor(TextbookOutcomeOR$DIALYSIS)
TextbookOutcomeOR$DISCANCR <- as.factor(TextbookOutcomeOR$DISCANCR)
TextbookOutcomeOR$HXCHF <- as.factor(TextbookOutcomeOR$HXCHF)
TextbookOutcomeOR$Age <- as.factor(TextbookOutcomeOR$Age)


TextbookOutcomeOR$SEX = relevel(TextbookOutcomeOR$SEX, ref="male")
TextbookOutcomeOR$BMILabled = relevel(TextbookOutcomeOR$BMILabled, ref="Normal Weight")
TextbookOutcomeOR$ASACLAS = relevel(TextbookOutcomeOR$ASACLAS, ref="4-Life Threat")
TextbookOutcomeOR$DIABETES = relevel(TextbookOutcomeOR$DIABETES, ref="Yes")
TextbookOutcomeOR$SMOKE = relevel(TextbookOutcomeOR$SMOKE, ref="Yes")
TextbookOutcomeOR$STEROID = relevel(TextbookOutcomeOR$STEROID, ref="Yes")
TextbookOutcomeOR$ASCITES = relevel(TextbookOutcomeOR$ASCITES, ref="Yes")
TextbookOutcomeOR$HXCOPD = relevel(TextbookOutcomeOR$HXCOPD, ref="Yes")
TextbookOutcomeOR$DIALYSIS = relevel(TextbookOutcomeOR$DIALYSIS, ref="Yes")
TextbookOutcomeOR$DISCANCR = relevel(TextbookOutcomeOR$DISCANCR, ref="Yes")
TextbookOutcomeOR$HXCHF = relevel(TextbookOutcomeOR$HXCHF, ref="Yes")
TextbookOutcomeOR$Age = relevel(TextbookOutcomeOR$Age, ref="85+")
####Multivariant analysis serious complications
MVORAll <- 
  glm(TextbookOutcomeAllOR ~ PatientGroup + Age + SEX + RACE_NEW + BMILabled + ASACLAS + STEROID + ASCITES + DYSPNEA + SMOKE + DIALYSIS + DISCANCR + DIABETES + HYPERMED + HXCOPD + HXCHF, 
      data = TextbookOutcomeOR,
      family = binomial("logit"), 
      na.action =na.omit)
###continue on to merge the table sets
MVORTableAll <-
  tbl_regression(MVORAll, exponentiate = T,
                 pvalue_fun = ~style_pvalue(.x, digits = 2),
  ) %>%
  bold_p(t = 0.10) %>%
  bold_labels() %>%
  italicize_levels() 
tbl_merge(
  list(UVAll, MVORTableAll),
  tab_spanner = c("**Univariable**", "**Multivariable**")
) %>%
  bold_labels() %>%
  italicize_levels()
Characteristic Univariable Multivariable
N OR1 95% CI1 p-value OR1 95% CI1 p-value
PatientGroup 35,154
Liver — — — —
Pancreas 0.77 0.74, 0.81 <0.001 0.72 0.68, 0.76 <0.001
Age 35,905
85+ — — — —
0-65 1.76 1.47, 2.12 <0.001 1.63 1.33, 2.00 <0.001
65-74 1.42 1.18, 1.72 <0.001 1.43 1.17, 1.76 <0.001
75-84 1.19 0.98, 1.44 0.075 1.23 1.00, 1.51 0.055
SEX 35,905
male — — — —
female 1.28 1.23, 1.33 <0.001 1.22 1.17, 1.28 <0.001
RACE_NEW 30,147
Black — — — —
Other 0.93 0.83, 1.04 0.20 0.81 0.72, 0.91 <0.001
White 1.11 1.03, 1.19 0.009 1.11 1.02, 1.20 0.012
BMILabled 35,723
Normal Weight — — — —
Morbidly Obese 0.79 0.71, 0.87 <0.001 0.74 0.66, 0.84 <0.001
Obese 0.87 0.82, 0.92 <0.001 0.86 0.81, 0.92 <0.001
Overweight 0.95 0.90, 1.00 0.056 0.96 0.91, 1.02 0.23
Underweight 0.96 0.83, 1.12 0.61 0.98 0.83, 1.16 0.83
ASACLAS 35,830
4-Life Threat — — — —
1-No Disturb 2.97 2.37, 3.73 <0.001 2.02 1.48, 2.78 <0.001
2-Mild Disturb 2.33 2.13, 2.56 <0.001 1.83 1.62, 2.07 <0.001
3-Severe Disturb 1.65 1.52, 1.81 <0.001 1.45 1.30, 1.63 <0.001
STEROID 35,905
Yes — — — —
No 1.31 1.16, 1.49 <0.001 1.21 1.06, 1.39 0.006
ASCITES 35,905
Yes — — — —
No 2.45 1.73, 3.55 <0.001 2.35 1.61, 3.52 <0.001
DYSPNEA 35,905
AT REST — — — —
MODERATE EXERTION 1.20 0.71, 2.09 0.52 1.14 0.64, 2.09 0.66
No 1.79 1.07, 3.11 0.032 1.42 0.81, 2.59 0.24
SMOKE 35,905
Yes — — — —
No 1.03 0.97, 1.09 0.36 1.02 0.95, 1.09 0.61
DIALYSIS 35,905
Yes — — — —
No 2.71 1.84, 4.12 <0.001 2.08 1.34, 3.35 0.002
DISCANCR 35,905
Yes — — — —
No 1.05 0.99, 1.10 0.091 1.24 1.16, 1.33 <0.001
DIABETES 35,905
Yes — — — —
NO 1.24 1.18, 1.30 <0.001 1.02 0.96, 1.09 0.44
HYPERMED 35,905
No — — — —
Yes 0.78 0.75, 0.81 <0.001 0.93 0.88, 0.98 0.005
HXCOPD 35,905
Yes — — — —
No 1.71 1.53, 1.91 <0.001 1.41 1.25, 1.60 <0.001
HXCHF 35,905
Yes — — — —
No 2.73 1.89, 4.04 <0.001 1.74 1.16, 2.69 0.010

1 OR = Odds Ratio, CI = Confidence Interval 

PanOddsPlotty <- 
  PanTextbookOutcomesSubset %>%
  select(PRNCPTX, SEX, RACE_NEW, BMILabled, STEROID, DIALYSIS, HXCOPD, HXCHF, TextbookOutcomePanOR)
PanOddsPlotty2 <- PanOddsPlotty[complete.cases(PanOddsPlotty), ] #Create a copy
head(PanOddsPlotty2, 10)
##                     PRNCPTX    SEX RACE_NEW     BMILabled STEROID DIALYSIS
## 17458               Whipple   male    White         Obese     Yes       No
## 17459               Whipple   male    White Normal Weight      No       No
## 17460               Whipple female    White         Obese      No       No
## 17461               Whipple female    White         Obese      No       No
## 17462               Whipple   male    White    Overweight      No       No
## 17463               Whipple female    Black         Obese      No       No
## 17464 Distal Pancreatectomy female    White Normal Weight     Yes       No
## 17465 Distal Pancreatectomy   male    White         Obese     Yes       No
## 17466 Distal Pancreatectomy   male    White         Obese      No       No
## 17467 Distal Pancreatectomy female    White         Obese      No       No
##       HXCOPD HXCHF TextbookOutcomePanOR
## 17458     No    No                    0
## 17459     No    No                    1
## 17460     No    No                    0
## 17461     No    No                    0
## 17462     No    No                    0
## 17463     No    No                    0
## 17464     No    No                    0
## 17465     No    No                    1
## 17466     No    No                    1
## 17467     No    No                    0
PanOddsPlotty2$TextbookOutcomePanOR <- factor(PanOddsPlotty2$TextbookOutcomePanOR)
PanOddsPlotty2$RACE_NEW <- factor(PanOddsPlotty2$RACE_NEW)

str(PanOddsPlotty2)
## 'data.frame':    19072 obs. of  9 variables:
##  $ PRNCPTX             : Factor w/ 3 levels "Distal Pancreatectomy",..: 3 3 3 3 3 3 1 1 1 1 ...
##  $ SEX                 : Factor w/ 2 levels "male","female": 1 1 2 2 1 2 2 1 1 2 ...
##  $ RACE_NEW            : Factor w/ 3 levels "Black","Other",..: 3 3 3 3 3 1 3 3 3 3 ...
##  $ BMILabled           : Factor w/ 5 levels "Normal Weight",..: 4 1 4 4 2 4 1 4 4 4 ...
##  $ STEROID             : Factor w/ 2 levels "Yes","No": 1 2 2 2 2 2 1 1 2 2 ...
##  $ DIALYSIS            : Factor w/ 2 levels "Yes","No": 2 2 2 2 2 2 2 2 2 2 ...
##  $ HXCOPD              : Factor w/ 2 levels "Yes","No": 2 2 2 2 2 2 2 2 2 2 ...
##  $ HXCHF               : Factor w/ 2 levels "Yes","No": 2 2 2 2 2 2 2 2 2 2 ...
##  $ TextbookOutcomePanOR: Factor w/ 2 levels "0","1": 1 2 1 1 1 1 1 2 2 1 ...
for(i in 1:8) {
  PanOddsPlotty2[, i] <- as.numeric(as.factor(PanOddsPlotty2[, i]))
}
###renaming variables
PanOddsPlotty2 <- rename(PanOddsPlotty2, 'No Dialysis' = DIALYSIS)
PanOddsPlotty2 <- rename(PanOddsPlotty2, 'Female' = SEX)
PanOddsPlotty2 <- rename(PanOddsPlotty2, 'White' = RACE_NEW)
PanOddsPlotty2 <- rename(PanOddsPlotty2, 'No Steroid' = STEROID)
PanOddsPlotty2 <- rename(PanOddsPlotty2, 'No COPD' = HXCOPD)
PanOddsPlotty2 <- rename(PanOddsPlotty2, 'No CHF' = HXCHF)
PanOddsPlotty2 <- rename(PanOddsPlotty2, 'Total Pancreatectomy' = PRNCPTX)
PanOddsPlotty2 <- rename(PanOddsPlotty2, 'Underweight' = BMILabled)
library(caret)
glm_model <- train(TextbookOutcomePanOR ~ .,
                   data = PanOddsPlotty2,
                   method = "glm",
                   family = "binomial")

summary(glm_model)
## 
## Call:
## NULL
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.1376  -1.0364  -0.9599   1.3069   2.0774  
## 
## Coefficients:
##                              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                  -5.40283    0.81631  -6.619 3.63e-11 ***
## `\\`Total Pancreatectomy\\``  0.02087    0.01599   1.305  0.19189    
## Female                        0.20488    0.02975   6.887 5.71e-12 ***
## White                         0.07510    0.02421   3.102  0.00192 ** 
## Underweight                  -0.07883    0.01199  -6.572 4.96e-11 ***
## `\\`No Steroid\\``            0.24694    0.08968   2.753  0.00590 ** 
## `\\`No Dialysis\\``           0.88736    0.29118   3.047  0.00231 ** 
## `\\`No COPD\\``               0.47938    0.07734   6.198 5.71e-10 ***
## `\\`No CHF\\``                0.73109    0.26413   2.768  0.00564 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 25810  on 19071  degrees of freedom
## Residual deviance: 25635  on 19063  degrees of freedom
## AIC: 25653
## 
## Number of Fisher Scoring iterations: 4
plotty <- OddsPlotty::odds_plot(glm_model$finalModel, 
                      title = "Pancreatic Surgery Group",
                      subtitle = "Independent Associations for TO")
## Waiting for profiling to be done...
plot <- plotty$odds_plot

plot <- plot + ggthemes::theme_economist() + theme(legend.position = "NULL")

plot + geom_text(label=round(plotty$odds_plot$data$OR, digits = 2), hjust=0, vjust=1.8)

HepOddsPlotty <- 
  HepTextbookOutcomesSubset %>%
  select(PRNCPTX, SEX, RACE_NEW, BMILabled, STEROID, DIALYSIS, DIABETES, HXCOPD, HXCHF, TextbookOutcomeHepOR)
HepOddsPlotty2 <- HepOddsPlotty[complete.cases(HepOddsPlotty), ] #Create a copy
head(HepOddsPlotty2, 10)
##              PRNCPTX    SEX RACE_NEW     BMILabled STEROID DIALYSIS DIABETES
## 1  Partial Lobectomy   male    White    Overweight      No       No      Yes
## 2  Partial Lobectomy female    White         Obese      No       No       NO
## 3  Partial Lobectomy female    Black         Obese      No       No       NO
## 4  Partial Lobectomy female    White Normal Weight      No       No       NO
## 5  Partial Lobectomy   male    White    Overweight      No       No       NO
## 6  Partial Lobectomy   male    White    Overweight      No       No       NO
## 7  Partial Lobectomy female    White         Obese      No       No       NO
## 8  Partial Lobectomy female    White         Obese      No       No       NO
## 9  Partial Lobectomy   male    White    Overweight      No       No      Yes
## 10 Partial Lobectomy   male    White    Overweight      No       No       NO
##    HXCOPD HXCHF TextbookOutcomeHepOR
## 1      No    No                    0
## 2      No    No                    1
## 3      No    No                    1
## 4      No    No                    1
## 5      No    No                    1
## 6     Yes    No                    0
## 7      No    No                    1
## 8      No    No                    1
## 9      No    No                    0
## 10     No    No                    1
HepOddsPlotty2$TextbookOutcomeHepOR <- factor(HepOddsPlotty2$TextbookOutcomeHepOR)
HepOddsPlotty2$RACE_NEW <- factor(HepOddsPlotty2$RACE_NEW)

str(HepOddsPlotty2)
## 'data.frame':    11617 obs. of  10 variables:
##  $ PRNCPTX             : Factor w/ 3 levels "Trisegmentectomy",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ SEX                 : Factor w/ 2 levels "male","female": 1 2 2 2 1 1 2 2 1 1 ...
##  $ RACE_NEW            : Factor w/ 3 levels "Black","Other",..: 3 3 1 3 3 3 3 3 3 3 ...
##  $ BMILabled           : Factor w/ 5 levels "Normal Weight",..: 4 3 3 1 4 4 3 3 4 4 ...
##  $ STEROID             : Factor w/ 2 levels "Yes","No": 2 2 2 2 2 2 2 2 2 2 ...
##  $ DIALYSIS            : Factor w/ 2 levels "Yes","No": 2 2 2 2 2 2 2 2 2 2 ...
##  $ DIABETES            : Factor w/ 2 levels "Yes","NO": 1 2 2 2 2 2 2 2 1 2 ...
##  $ HXCOPD              : Factor w/ 2 levels "Yes","No": 2 2 2 2 2 1 2 2 2 2 ...
##  $ HXCHF               : Factor w/ 2 levels "Yes","No": 2 2 2 2 2 2 2 2 2 2 ...
##  $ TextbookOutcomeHepOR: Factor w/ 2 levels "0","1": 1 2 2 2 2 1 2 2 1 2 ...
for(i in 1:9) {
  HepOddsPlotty2[, i] <- as.numeric(as.factor(HepOddsPlotty2[, i]))
}
###renaming variables
HepOddsPlotty2 <- rename(HepOddsPlotty2, 'No Dialysis' = DIALYSIS)
HepOddsPlotty2 <- rename(HepOddsPlotty2, 'Female' = SEX)
HepOddsPlotty2 <- rename(HepOddsPlotty2, 'White' = RACE_NEW)
HepOddsPlotty2 <- rename(HepOddsPlotty2, 'No Steroid' = STEROID)
HepOddsPlotty2 <- rename(HepOddsPlotty2, 'No COPD' = HXCOPD)
HepOddsPlotty2 <- rename(HepOddsPlotty2, 'No CHF' = HXCHF)
HepOddsPlotty2 <- rename(HepOddsPlotty2, 'No Diabetes' = DIABETES)

HepOddsPlotty2 <- rename(HepOddsPlotty2, 'Lobectomy' = PRNCPTX)
HepOddsPlotty2 <- rename(HepOddsPlotty2, 'Obese' = BMILabled)
library(caret)
glm_modelhep <- train(TextbookOutcomeHepOR ~ .,
                   data = HepOddsPlotty2,
                   method = "glm",
                   family = "binomial")

summary(glm_modelhep)
## 
## Call:
## NULL
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.3185  -1.1638  -0.8909   1.1726   1.9965  
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -6.49386    1.04277  -6.227 4.74e-10 ***
## Lobectomy            0.11848    0.03004   3.944 8.01e-05 ***
## Female               0.28058    0.03775   7.432 1.07e-13 ***
## White                0.04367    0.02887   1.512  0.13043    
## Obese                0.03377    0.01508   2.240  0.02512 *  
## `\\`No Steroid\\``   0.19966    0.10526   1.897  0.05786 .  
## `\\`No Dialysis\\``  0.83174    0.37348   2.227  0.02595 *  
## `\\`No Diabetes\\``  0.33922    0.04961   6.838 8.03e-12 ***
## `\\`No COPD\\``      0.43853    0.09746   4.500 6.81e-06 ***
## `\\`No CHF\\``       0.99246    0.34857   2.847  0.00441 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 16100  on 11616  degrees of freedom
## Residual deviance: 15918  on 11607  degrees of freedom
## AIC: 15938
## 
## Number of Fisher Scoring iterations: 4
plottyhep <- OddsPlotty::odds_plot(glm_modelhep$finalModel, 
                      title = "Hepatic Surgery Group",
                      subtitle = "Independent Associations for TO")
## Waiting for profiling to be done...
plothep <- plottyhep$odds_plot

plothep <- plothep + ggthemes::theme_economist() + theme(legend.position = "NULL")

plothep + geom_text(label=round(plottyhep$odds_plot$data$OR, digits = 2), hjust=0, vjust=1.8)

BothOddsPlotty <- 
  TextbookOutcomeOR %>%
  select(PatientGroup, SEX, STEROID, DIALYSIS, DIABETES, HXCOPD, HXCHF, TextbookOutcomeAllOR)
BothOddsPlotty2 <- BothOddsPlotty[complete.cases(BothOddsPlotty), ] #Create a copy
head(BothOddsPlotty2, 10)
##    PatientGroup    SEX STEROID DIALYSIS DIABETES HXCOPD HXCHF
## 1         Liver   male      No       No      Yes     No    No
## 2         Liver female      No       No       NO     No    No
## 3         Liver female      No       No       NO     No    No
## 4         Liver female      No       No       NO     No    No
## 5         Liver   male      No       No       NO     No    No
## 6         Liver   male      No       No       NO    Yes    No
## 7         Liver female      No       No       NO     No    No
## 8         Liver female      No       No       NO     No    No
## 9         Liver   male      No       No      Yes     No    No
## 10        Liver   male      No       No       NO     No    No
##    TextbookOutcomeAllOR
## 1                     0
## 2                     1
## 3                     1
## 4                     1
## 5                     1
## 6                     0
## 7                     1
## 8                     1
## 9                     0
## 10                    1
BothOddsPlotty2$TextbookOutcomeAllOR <- factor(BothOddsPlotty2$TextbookOutcomeAllOR)

str(BothOddsPlotty2)
## 'data.frame':    35154 obs. of  8 variables:
##  $ PatientGroup        : Factor w/ 2 levels "Liver","Pancreas": 1 1 1 1 1 1 1 1 1 1 ...
##  $ SEX                 : Factor w/ 2 levels "male","female": 1 2 2 2 1 1 2 2 1 1 ...
##  $ STEROID             : Factor w/ 2 levels "Yes","No": 2 2 2 2 2 2 2 2 2 2 ...
##  $ DIALYSIS            : Factor w/ 2 levels "Yes","No": 2 2 2 2 2 2 2 2 2 2 ...
##  $ DIABETES            : Factor w/ 2 levels "Yes","NO": 1 2 2 2 2 2 2 2 1 2 ...
##  $ HXCOPD              : Factor w/ 2 levels "Yes","No": 2 2 2 2 2 1 2 2 2 2 ...
##  $ HXCHF               : Factor w/ 2 levels "Yes","No": 2 2 2 2 2 2 2 2 2 2 ...
##  $ TextbookOutcomeAllOR: Factor w/ 2 levels "0","1": 1 2 2 2 2 1 2 2 1 2 ...
for(i in 1:7) {
  BothOddsPlotty2[, i] <- as.numeric(as.factor(BothOddsPlotty2[, i]))
}
###renaming variables
BothOddsPlotty2 <- rename(BothOddsPlotty2, 'No Dialysis' = DIALYSIS)
BothOddsPlotty2 <- rename(BothOddsPlotty2, 'Female' = SEX)
BothOddsPlotty2 <- rename(BothOddsPlotty2, 'No Steroid' = STEROID)
BothOddsPlotty2 <- rename(BothOddsPlotty2, 'No COPD' = HXCOPD)
BothOddsPlotty2 <- rename(BothOddsPlotty2, 'No CHF' = HXCHF)
BothOddsPlotty2 <- rename(BothOddsPlotty2, 'No Diabetes' = DIABETES)
BothOddsPlotty2 <- rename(BothOddsPlotty2, 'Pancreatic' = PatientGroup)
library(caret)
glm_modelboth <- train(TextbookOutcomeAllOR ~ .,
                   data = BothOddsPlotty2,
                   method = "glm",
                   family = "binomial")

summary(glm_modelboth)
## 
## Call:
## NULL
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.2248  -1.1220  -0.9696   1.2338   2.0334  
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -5.36928    0.59904  -8.963  < 2e-16 ***
## Pancreatic          -0.24256    0.02201 -11.022  < 2e-16 ***
## Female               0.22728    0.02167  10.488  < 2e-16 ***
## `\\`No Steroid\\``   0.24218    0.06422   3.771 0.000162 ***
## `\\`No Dialysis\\``  0.90930    0.21287   4.272 1.94e-05 ***
## `\\`No Diabetes\\``  0.15174    0.02612   5.809 6.27e-09 ***
## `\\`No COPD\\``      0.48508    0.05738   8.453  < 2e-16 ***
## `\\`No CHF\\``       0.84574    0.19760   4.280 1.87e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 48336  on 35153  degrees of freedom
## Residual deviance: 47891  on 35146  degrees of freedom
## AIC: 47907
## 
## Number of Fisher Scoring iterations: 4
plottyboth <- OddsPlotty::odds_plot(glm_modelboth$finalModel, 
                      title = "Combined Liver & Pancreatic Groups",
                      subtitle = "Independent Associations for TO")
## Waiting for profiling to be done...
plotboth <- plottyboth$odds_plot

plotboth <- plotboth + ggthemes::theme_economist() + theme(legend.position = "NULL")

plotboth + geom_text(label=round(plottyboth$odds_plot$data$OR, digits = 2), hjust=0, vjust=1.8)