Adrenalectomy

Ibrahim Hassan

2023-09-20

library(gtsummary)
library(tidyverse)
library(readxl)
library(skimr)
library(ggpubr)
library(flextable)
library(gt)
library(scales)
library(xtable)
library(pillar)
library(Hmisc)
library(kableExtra)
library(paletteer)
library(ggthemes)
library(ggsci)
library(hrbrthemes)
library(easystats)
library(reshape2)
df <- read_excel("C:/Users/HiKa8/Desktop/Work/Adrenalectomy/Adrenelectomy_data_Sept29_clean.xlsx" )
df$Clavien <- factor(df$`Clavien-Dindo`, levels = c("0", "1", "2", "3", "4", "5"))
df$Surgery_Service <- factor(df$Surgery_Service)
df$Indication <- factor(df$Indication)
df$`30day_readmit` <- factor(df$`30day_readmit`)
df$ICU_admission <- factor(df$ICU_admission, levels = c("Yes/Step-down ICU", "No"))

Patient statistics

Patient characteristics

df %>%
  select(Age,Surgery_Service, Charlson_Comorbidity ) %>%
  tbl_summary(
    by = Surgery_Service,
    type = all_continuous() ~ "continuous2",
    digits = list(all_continuous2() ~ 2, 
                  all_categorical() ~ c(0,2)),
    statistic = all_continuous2() ~ c("{min} - {max}",
                                     "{mean} (±{sd})",
                                     "{median} ({p25}, {p75})"),
    label = list(Charlson_Comorbidity ~ "Charlson Comorbidity score")
  ) %>%
  add_p(list(
      all_continuous2() ~ "t.test",
      all_categorical() ~ "chisq.test"
      )
  ) %>%
  bold_p( t = 0.05, q = FALSE)
Characteristic GENERAL, N = 37 UROLOGY, N = 66 p-value1
Age 0.3
    Range 13.00 - 78.00 26.00 - 79.00
    Mean (±SD) 53.92 (±16.24) 57.17 (±13.52)
    Median (IQR) 56.00 (47.00, 67.00) 60.00 (50.25, 66.75)
Charlson Comorbidity score 0.7
    Range 0.00 - 10.00 2.00 - 12.00
    Mean (±SD) 4.62 (±2.58) 4.80 (±2.47)
    Median (IQR) 4.00 (3.00, 6.00) 4.50 (3.00, 6.00)
1 Welch Two Sample t-test

Histogram showing the age disribution categorized by surgery service

gghistogram(df, x = "Age",
   add = "mean", rug = TRUE,
   fill = "Surgery_Service",
   add_density = TRUE) +
  scale_fill_wsj()

Difference of patient characteristics by surgery type

df %>% 
   ggboxplot(
    y = "Charlson_Comorbidity", 
    rug = TRUE, 
    x = "Surgery_Service", 
    fill = "Surgery_Service", 
    ylab = "Charlson Comorbidity score",
    xlab = "",
    title = "", 
    palette = "npg", 
    ggtheme = theme_pubr(),bxp.errorbar = TRUE
  ) +
  stat_compare_means(method = "t.test")+
  rremove("legend") + 
  scale_fill_jama()+
  labs_pubr()

Surgery Statistics

Surgery characteristics

df %>%
  select(Surgery_Service,minutes, Indication, minutes ,`Open/Laparoscopic` ) %>%
  tbl_summary(
    by = Surgery_Service,
    type = all_continuous() ~ "continuous2",
    digits = list(all_continuous2() ~ 1, 
                  all_categorical() ~ c(0,2)),
    statistic = all_continuous2() ~ c("{min} - {max}",
                                     "{mean} (±{sd})",
                                     "{median} ({p25}, {p75})"),
    label = list(minutes ~ "Operating Room time (in minutes)"),
    sort = list(Indication ~ "frequency")
  ) %>%
  add_p(list(
      all_continuous2() ~ "t.test",
      all_categorical() ~ "chisq.test"
      )
  ) %>%
  bold_p( t = 0.05, q = FALSE)
Characteristic GENERAL, N = 371 UROLOGY, N = 661 p-value2
Operating Room time (in minutes) 0.7
    Range 158.0 - 618.0 123.0 - 490.0
    Mean (±SD) 271.8 (±111.5) 263.1 (±75.3)
    Median (IQR) 235.0 (206.0, 296.0) 256.0 (211.0, 302.3)
Indication 0.2
    Pheochromocytoma 12 (32.43%) 20 (30.30%)
    Functional Adrenal Lesion 5 (13.51%) 19 (28.79%)
    Metastases 9 (24.32%) 10 (15.15%)
    Size/growth 5 (13.51%) 11 (16.67%)
    Primary Adrenocortical Carcinoma 4 (10.81%) 6 (9.09%)
    Other 2 (5.41%) 0 (0.00%)
Open/Laparoscopic 0.4
    Laparoscopic 31 (83.78%) 61 (92.42%)
    Open 5 (13.51%) 4 (6.06%)
    Open/Laparoscopic 1 (2.70%) 1 (1.52%)
1 n (%)
2 Welch Two Sample t-test; Pearson’s Chi-squared test

Association of Operating Room time (in minutes) with Surgery service

df %>% 
   ggboxplot(
    y = "minutes",
    rug = TRUE, 
    x = "Surgery_Service", 
    fill = "Surgery_Service", 
    ylab = "Operating Room time (in minutes)",
    xlab = "Surgery Service",
    title = "", 
    palette = "npg", 
    ggtheme = theme_pubr(),
    bxp.errorbar = TRUE
  ) +
  stat_compare_means(label.x.npc = "middle", method = "t.test")+
  rremove("legend") +
  scale_fill_jama()+
  labs_pubr()+
  theme(plot.title = element_text(hjust = .9, vjust=.12))

Association of Operating Room time (in minutes) with Indication

df %>% 
   ggboxplot(
    y = "minutes",
    rug = TRUE, 
    x = "Indication", 
    fill = "Indication", 
    ylab = "Operating Room time (in minutes)",
    xlab = "",
    title = "", 
    ggtheme = theme_pubr(),
    bxp.errorbar = TRUE,
    orientation = "horizontal"
  ) +
  stat_compare_means( method = "anova", label.y.npc = .7)+
  rremove("legend") +
  labs_pubr()+
  theme(plot.title = element_text(hjust = .9, vjust=.12))

Surgical outcomes statistics

Surgical outcomes characterstics

df %>%
  select(Surgery_Service,`30day_readmit`, ICU_admission, Clavien, `Postoperative stay (days)`) %>%
  tbl_summary(
    by = Surgery_Service,
    type = all_continuous() ~ "continuous2",
    digits = list(all_continuous2() ~ 2, 
                  all_categorical() ~ c(0,2)),
    statistic = all_continuous2() ~ c("{min} - {max}",
                                     "{mean} (±{sd})",
                                     "{median} ({p25}, {p75})"),
   label = list(
     Clavien ~ "Clavien-Dindo score",
     `30day_readmit` ~ "30-day readmission",
     ICU_admission ~ "ICU admission"
     )  
  ) %>%
  add_p(
    list(
      all_continuous() ~ "t.test",
      all_categorical() ~ "chisq.test"
      )
    ) %>%
  bold_p( t = 0.05, q = FALSE)
Characteristic GENERAL, N = 371 UROLOGY, N = 661 p-value2
30-day readmission 5 (13.51%) 1 (1.52%) 0.040
ICU admission 0.039
    Yes/Step-down ICU 19 (51.35%) 19 (28.79%)
    No 18 (48.65%) 47 (71.21%)
Clavien-Dindo score 0.4
    0 23 (62.16%) 49 (74.24%)
    1 6 (16.22%) 8 (12.12%)
    2 5 (13.51%) 6 (9.09%)
    3 2 (5.41%) 1 (1.52%)
    4 0 (0.00%) 2 (3.03%)
    5 1 (2.70%) 0 (0.00%)
Postoperative stay (days) 0.068
    Range 1.00 - 50.00 1.00 - 9.00
    Mean (±SD) 5.00 (±8.05) 2.48 (±1.64)
    Median (IQR) 3.00 (2.00, 4.00) 2.00 (1.00, 3.00)
1 n (%)
2 Pearson’s Chi-squared test; Welch Two Sample t-test

Association of Surgery survice with 30 day readmission

chisq.test(df$Surgery_Service, df$`30day_readmit`) 
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  df$Surgery_Service and df$`30day_readmit`
## X-squared = 4.2267, df = 1, p-value = 0.03979
df %>%
  select(Surgery_Service,`30day_readmit`) %>%
  tbl_summary(
    by = Surgery_Service,
    type = all_continuous() ~ "continuous2",
    digits = list(all_continuous2() ~ 2, 
                  all_categorical() ~ c(0,2)),
    statistic = all_continuous2() ~ c("{min} - {max}",
                                     "{mean} (±{sd})",
                                     "{median} ({p25}, {p75})"),
    
   label = list(
     `30day_readmit` ~ "30-day readmission"
     )  
  ) %>%
  add_p(
    list(
      all_continuous() ~ "t.test",
      all_categorical() ~ "chisq.test"
      )
    ) %>%
  bold_p( t = 0.05, q = FALSE)
Characteristic GENERAL, N = 371 UROLOGY, N = 661 p-value2
30-day readmission 5 (13.51%) 1 (1.52%) 0.040
1 n (%)
2 Pearson’s Chi-squared test

Association of Postoperative stay with surgery service

t.test(`Postoperative stay (days)` ~ Surgery_Service, data = df) %>% report()
## Effect sizes were labelled following Cohen's (1988) recommendations.
## 
## The Welch Two Sample t-test testing the difference of Postoperative stay (days)
## by Surgery_Service (mean in group GENERAL = 5.00, mean in group UROLOGY = 2.48)
## suggests that the effect is positive, statistically not significant, and medium
## (difference = 2.52, 95% CI [-0.19, 5.22], t(37.68) = 1.88, p = 0.068; Cohen's d
## = 0.61, 95% CI [-0.04, 1.26])
df %>%
  select(Surgery_Service, `Postoperative stay (days)`) %>%
  tbl_summary(
    by = Surgery_Service,
    type = all_continuous() ~ "continuous2",
    digits = list(all_continuous2() ~ 2, 
                  all_categorical() ~ c(0,2)),
    statistic = all_continuous2() ~ c("{min} - {max}",
                                     "{mean} (±{sd})",
                                     "{median} ({p25}, {p75})")
  ) %>%
  add_p(
    list(
      all_continuous() ~ "t.test",
      all_categorical() ~ "chisq.test"
      )
    ) %>%
  bold_p( t = 0.05, q = FALSE)
Characteristic GENERAL, N = 37 UROLOGY, N = 66 p-value1
Postoperative stay (days) 0.068
    Range 1.00 - 50.00 1.00 - 9.00
    Mean (±SD) 5.00 (±8.05) 2.48 (±1.64)
    Median (IQR) 3.00 (2.00, 4.00) 2.00 (1.00, 3.00)
1 Welch Two Sample t-test

Association of ICU admission with surgery service

chisq.test(df$Surgery_Service, df$ICU_admission)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  df$Surgery_Service and df$ICU_admission
## X-squared = 4.2605, df = 1, p-value = 0.03901
df %>%
  select(Surgery_Service,ICU_admission) %>%
  tbl_summary(
    by = Surgery_Service,
    type = all_continuous() ~ "continuous2",
    digits = list(all_continuous2() ~ 2, 
                  all_categorical() ~ c(0,2)),
    statistic = all_continuous2() ~ c("{min} - {max}",
                                     "{mean} (±{sd})",
                                     "{median} ({p25}, {p75})"),
    label = list(ICU_admission ~ "ICU admission")
  ) %>%
  add_p(list(all_categorical() ~ "chisq.test")) %>%
  bold_p( t = 0.05, q = FALSE)
Characteristic GENERAL, N = 371 UROLOGY, N = 661 p-value2
ICU admission 0.039
    Yes/Step-down ICU 19 (51.35%) 19 (28.79%)
    No 18 (48.65%) 47 (71.21%)
1 n (%)
2 Pearson’s Chi-squared test

Frequency of multidisciplinary involvement (Bar plot)

x1 <- melt(x, id = c("Multidisciplinaryi_involvement"))[-1]

x1 %>% 
  group_by(variable) %>% 
  summarise(count = round(sum(value))) %>% 
  ggplot(aes(x = reorder(variable, (count)), y = count, fill = variable)) + 
  geom_bar(stat = 'identity') +
  coord_flip() +
  xlab("") +
  ylab("") +
  theme_pubr(legend = "none") +
  labs_pubr() +
  theme(plot.title = element_text(hjust = 3, vjust=.12)) 

x1 %>%
  group_by(variable) %>% 
  summarise(count = sum(value)) %>%
  rename("Descipline" = "variable") %>% 
  arrange(desc(count)) %>% 
  kbl(fromat = "html")
Descipline count
Endocrinology 64
Oncology 19
Internal Medicine 11
General surgery 4
Thoracic surgery 4
Urology 2
Palliative care 1
Respirology 1
Nephrology 1
Cardiology 1
Cardiac surgery 1
Gynecological-Oncology 1
Dermatology 1
ICU 1

Association of Clavien-Dindo score with Charlson Comorbidity score

df %>% 
   ggboxplot(
    y = "Charlson_Comorbidity",
    rug = TRUE, 
    x = "Clavien", 
    fill = "Clavien", 
    ylab = "Charlson Comorbidity score",
    xlab = "Clavien-Dindo score",
    title = "", 
    ggtheme = theme_pubr(),bxp.errorbar = TRUE
    
  ) +
  stat_compare_means(label.x.npc = "middle", method = "anova")+
  rremove("legend") + 
  scale_fill_futurama()+
  labs_pubr()+
  theme(plot.title = element_text(hjust = .9, vjust=.12))

Association of Age with Clavien-Dindo score

df %>% 
   ggboxplot(
    y = "Age",
    rug = TRUE, # adds a rug to the plot
    x = "Clavien", 
    fill = "Clavien", 
    ylab = "Age",
    xlab = "Clavien-Dindo score",
    title = "", 
    palette = "npg", 
    ggtheme = theme_pubr(),
    bxp.errorbar = TRUE
  ) +
  stat_compare_means(label.x.npc = "middle", method = "anova")+
  rremove("legend") +
  scale_fill_jama()+
  labs_pubr()+
  theme(plot.title = element_text(hjust = .9, vjust=.12))

Difference in Laparoscopic Vs Open Open Operating Room time (in minutes) according to different surgery services

df %>% 
    select(Surgery_Service, minutes,`Open/Laparoscopic`) %>%
    filter(`Open/Laparoscopic` == "Laparoscopic" | `Open/Laparoscopic` == "Open"  ) %>% 
    tbl_strata(
        strata = `Open/Laparoscopic`,
        .tbl_fun =
            ~ .x %>%
            tbl_summary(
                by = Surgery_Service,
                type = list(minutes ~ "continuous2"), 
                statistic = list(                 
                    all_continuous() ~ c("{min} - {max}",
                                         "{mean} (±{sd})",
                                         "{median} ({p25}, {p75})"
                                         ), 
                    all_categorical() ~ "{n} ({p}%)"
                ),
                digits = list(all_categorical() ~ c(0, 1)),
                missing_text = "No Response",
                label = list(minutes = "Open Operating Room time (in minutes)")
            ) %>% 
          add_p(list(all_continuous2() ~ "t.test"))
    )
Characteristic Laparoscopic Open
GENERAL, N = 31 UROLOGY, N = 61 p-value1 GENERAL, N = 5 UROLOGY, N = 4 p-value1
Open Operating Room time (in minutes) >0.9 0.8
    Range 158 - 618 123 - 487 209 - 548 271 - 361
    Mean (±SD) 253 (±98) 256 (±70) 335 (±127) 318 (±39)
    Median (IQR) 229 (201, 274) 251 (209, 299) 306 (289, 325) 319 (295, 342)
1 Welch Two Sample t-test

Difference in Laparoscopic Vs Open Operating Room time (in minutes) according to Difference surgeons

df %>% 
    select(Surgery_Service, minutes, Indication) %>%
    filter(Indication %in% c("Functional Adrenal Lesion", "Pheochromocytoma", "Size/growth")) %>% 
    tbl_strata(
      strata = Indication,
      .tbl_fun = ~ .x %>%
        tbl_summary(
            by = Surgery_Service,
            type = list(all_continuous() ~ "continuous2"), 
            
            statistic = list(                 
                all_continuous() ~ c("{min} - {max}",
                                     "{mean} (±{sd})",
                                     "{median} ({p25}, {p75})"
                                     ), 
                all_categorical() ~ "{n} ({p}%)"
            ),
            digits = list(all_categorical() ~ c(0, 1)),
            missing_text = "No Response",
            label = list(minutes ~ "Operating Room time (in minutes)")
        ) %>% 
        add_p(list(
            all_continuous2() ~ "t.test",
            all_categorical() ~ "chisq.test"
            )
        ) %>%
        bold_p( t = 0.05, q = FALSE)
    )
Characteristic Functional Adrenal Lesion Pheochromocytoma Size/growth
GENERAL, N = 5 UROLOGY, N = 19 p-value1 GENERAL, N = 12 UROLOGY, N = 20 p-value1 GENERAL, N = 5 UROLOGY, N = 11 p-value1
Operating Room time (in minutes) 0.5 0.005 0.7
    Range 219 - 255 142 - 394 172 - 309 199 - 407 158 - 548 123 - 328
    Mean (±SD) 232 (±15) 243 (±63) 238 (±46) 292 (±52) 257 (±163) 229 (±67)
    Median (IQR) 229 (220, 235) 249 (216, 261) 234 (202, 284) 293 (258, 317) 190 (185, 206) 233 (186, 271)
1 Welch Two Sample t-test

Association of Age with 30 day readmission

t.test(Age ~ `30day_readmit`, data = df) %>% report()
## Effect sizes were labelled following Cohen's (1988) recommendations.
## 
## The Welch Two Sample t-test testing the difference of Age by 30day_readmit
## (mean in group No = 56.10, mean in group Yes = 54.33) suggests that the effect
## is positive, statistically not significant, and small (difference = 1.77, 95%
## CI [-10.91, 14.45], t(5.95) = 0.34, p = 0.744; Cohen's d = 0.28, 95% CI [-1.34,
## 1.88])
df %>% 
  select(`30day_readmit`, Age) %>% 
  tbl_summary(
    by = `30day_readmit`,
    type = list(all_continuous() ~ "continuous2"), 
    statistic = list(                 
                all_continuous() ~ c("{min} - {max}",
                                    "{mean} (±{sd})",
                                    "{median} ({p25}, {p75})"
                                         ), 
                    all_categorical() ~ "{n} ({p}%)"
                ),
                digits = list(all_categorical() ~ c(0, 1)),
                missing_text = "No Response"
            ) %>% 
  add_p(list(
      all_continuous2() ~ "t.test",
      all_categorical() ~ "chisq.test"
      )
  ) %>%
  bold_p( t = 0.05, q = FALSE)
Characteristic No, N = 97 Yes, N = 6 p-value1
Age 0.7
    Range 13 - 79 35 - 67
    Mean (±SD) 56 (±15) 54 (±12)
    Median (IQR) 59 (49, 67) 55 (50, 64)
1 Welch Two Sample t-test

Association of Operating Room Time (in minutes) with 30 day readmission

t.test(minutes ~ `30day_readmit`, data = df) %>% report()
## Effect sizes were labelled following Cohen's (1988) recommendations.
## 
## The Welch Two Sample t-test testing the difference of minutes by 30day_readmit
## (mean in group No = 260.51, mean in group Yes = 358.50) suggests that the
## effect is negative, statistically not significant, and large (difference =
## -97.99, 95% CI [-276.91, 80.92], t(5.14) = -1.40, p = 0.220; Cohen's d = -1.23,
## 95% CI [-3.06, 0.69])
df %>% 
  select(`30day_readmit`,minutes) %>% 
    tbl_summary(
          by = `30day_readmit`,
          type = list(all_continuous() ~ "continuous2"), 
          
          statistic = list(                 
              all_continuous() ~ c("{min} - {max}",
                                   "{mean} (±{sd})",
                                   "{median} ({p25}, {p75})"
                                   ), 
              all_categorical() ~ "{n} ({p}%)"
          ),
          label = list(minutes ~ "Operating Room time (in minutes)"),
          digits = list(all_categorical() ~ c(0, 1)),
          missing_text = "No Response"
            ) %>% 
  modify_header(label = "**30-day readmission**") %>% 
    add_p(list(
      all_continuous2() ~ "t.test",
      all_categorical() ~ "chisq.test"
      )
  ) %>%
  bold_p( t = 0.05, q = FALSE)
30-day readmission No, N = 97 Yes, N = 6 p-value1
Operating Room time (in minutes) 0.2
    Range 123 - 548 229 - 618
    Mean (±SD) 261 (±80) 359 (±171)
    Median (IQR) 251 (206, 300) 271 (236, 472)
1 Welch Two Sample t-test

Association of Operating Room time (in minutes) with indication

aov(minutes ~ Indication, data = df) %>% report()
## The ANOVA (formula: minutes ~ Indication) suggests that:
## 
##   - The main effect of Indication is statistically not significant and medium
## (F(5, 97) = 2.18, p = 0.062; Eta2 = 0.10, 95% CI [0.00, 1.00])
## 
## Effect sizes were labelled following Field's (2013) recommendations.