LEAN

Code 
knitr::opts_chunk$set(fig.width = 14)
library(pacman)
p_load(tidyverse,
       gtsummary,
       janitor)
ggplot2::theme_set(theme_minimal())

Outline

This document is organized as follows:

  1. Section 1: Data loading - the data are read (from a preprocessed excel file) and cleaned.

  2. Section 2: Basic statistics - the available variables (excluding durations) are analyzed and compared across surgeons.

  3. Section 3: Analysis of surgery times - the surgery time is analyzed across surgeons, phases and their interaction.

  4. Section 4: New chapter (07-11-22) dealing with each surgeon’s personal performance.

  5. ?@sec-hypotheses: Inference about the outcomes of phase-duration + “testing” (or cancelling) your (supervisor’s) original hypotheses + additional analyses.

  6. ?@sec-recommendations: Recommendations - for the current reasearch project and future works.

Data loading & cleaning

  • Zepto was “yes” only in one observation - variable removed.

  • Missing complexity was recoded to NA.

Code 
set.seed(1312)
lean <- read_csv("LEAN_final_210922.csv") %>% 
  janitor::clean_names() %>% 
  select(-surgery_date, -file_name, -zepto)
  

lean$complexity <- 
  parse_number(lean$complexity ) %>% as.vector()

lean$surgeon <- paste0("Surgeon",lean$surgeon)

Basic statistics

Descriptive statistics are provided for all study variables, overall and by surgeon. Remarks of interesting and important findings are given in the bullets below.

Phaco technique

  • Surgeon 1 “prefers”1 the partial divide-and-conquer technique (used in 35/50 of his surgeries).

  • Surgeon 2 has employed almost exlusively (47 of 50) the horizontal-chop technique.

  • For Surgeon 3, horizontal chop is the predominant technique (40 of 50).

  • Surgeon 4 uses the horizontal chop and the stop-n-chop techniques in most of her surgeries (42% and 38%, respectively).

  • Surgeon 1 is a non-conformist who popularizes the uncommon partial divide-and-conquer technique - note how 7 out of 8 surgeries performed in that manner were done by surgeon 1.

Keep the following questions in mind:

Code 
lean %>% 
  select(phaco_technique,surgeon) %>% 
  tbl_summary(by = surgeon) %>% 
  add_overall()
Characteristic Overall, N = 2001 Surgeon1, N = 501 Surgeon2, N = 501 Surgeon3, N = 501 Surgeon4, N = 501
phaco_technique
    divide and conquer 3 (1.5%) 1 (2.0%) 0 (0%) 0 (0%) 2 (4.0%)
    horizontal chop 110 (55%) 2 (4.0%) 47 (94%) 40 (80%) 21 (42%)
    partial divide and conquer 40 (20%) 35 (70%) 0 (0%) 0 (0%) 5 (10%)
    pre chopper 4 (2.0%) 0 (0%) 0 (0%) 4 (8.0%) 0 (0%)
    stop and chop 33 (17%) 11 (22%) 1 (2.0%) 2 (4.0%) 19 (38%)
    vertical chop 10 (5.0%) 1 (2.0%) 2 (4.0%) 4 (8.0%) 3 (6.0%)
1 n (%)

IOL Injection

  • Evidently, surgeons 2 and 3 focused on “BSS” injection almost exclusively, while surgeon 4 performed 10 (out of 48) “OVD” injections. In contrast, the avant-garde surgeon 1 performed 20 (40%) OVD injections.
Code 
lean %>% 
  select(iol_injection_under,surgeon) %>% 
  tbl_summary(by = surgeon) %>% 
  add_overall()
Characteristic Overall, N = 2001 Surgeon1, N = 501 Surgeon2, N = 501 Surgeon3, N = 501 Surgeon4, N = 501
iol_injection_under
    BSS 167 (84%) 30 (60%) 48 (96%) 50 (100%) 39 (78%)
    OVD 33 (17%) 20 (40%) 2 (4.0%) 0 (0%) 11 (22%)
1 n (%)

Duration statistics

In this section we analyze the duration, across operation phases and surgeons.

By phase

Descriptive statistics for each operational phase:

Code 
dur <- 
lean %>% 
  select(ends_with("duration")) %>% 
  rename_with(~str_remove(.x,"_duration")) %>% 
  rename_with(~str_replace_all(.x,"_"," "))


times_table <- dur %>% tbl_summary(
                  type = all_continuous() ~ "continuous2",
    statistic = all_continuous() ~ c("{mean} ({sd})", "{median} ({p25}, {p75})", "{min}, {max}")) 

times_table
Characteristic N = 200
surgical ports
    Mean (SD) 97 (44)
    Median (IQR) 91 (72, 111)
    Range 34, 271
capsulorhexis
    Mean (SD) 74 (32)
    Median (IQR) 67 (55, 84)
    Range 33, 336
hydrodissection and hydrodelineation
    Mean (SD) 71 (46)
    Median (IQR) 59 (42, 85)
    Range 23, 376
phacoemulsification
    Mean (SD) 282 (141)
    Median (IQR) 243 (191, 334)
    Range 87, 1,072
cortex aspiration
    Mean (SD) 140 (77)
    Median (IQR) 125 (84, 173)
    Range 43, 578
iol preparation
    Mean (SD) 92 (51)
    Median (IQR) 87 (64, 109)
    Range 13, 376
iol implantation
    Mean (SD) 70 (54)
    Median (IQR) 55 (40, 80)
    Range 15, 450
hydrating the ports and injection of antibiotic
    Mean (SD) 78 (70)
    Median (IQR) 57 (34, 93)
    Range 9, 521
total surgery
    Mean (SD) 905 (342)
    Median (IQR) 829 (659, 1,055)
    Range 432, 2,268

Histograms of phase duration:

Code 
durlong <- dur %>% 
  pivot_longer(everything()) %>% 
  rename(phase = name, duration = value) %>% 
  mutate(phase = fct(phase))


durlong$phase <- 
  fct_recode(durlong$phase,
             "Surgical ports" = "surgical ports",
             "Capsulorhexis" = "capsulorhexis",
             "Hydrodissection & hydrodelineation" = "hydrodissection and hydrodelineation",
             "Phacoemulsification" = "phacoemulsification",
             "Cortex aspiration" = "cortex aspiration",
             "IOL preparation" = "iol preparation",
             "IOL implantation" = "iol implantation",
             "Hydrating the ports & ABX injection" = "hydrating the ports and injection of antibiotic",
             "Total surgery" = "total surgery")
Code 
durlong %>% 
  ggplot(aes(x=duration))+
  facet_wrap(~phase,scales = "free")+geom_histogram()+
  jtools::theme_apa()+
  ggtitle("Duration by phase","units = [seconds]")

Same histograms, with duration log transformed:

Code 
durlong %>% 
  ggplot(aes(x=log(duration)))+
  facet_wrap(~phase,scales = "free")+geom_histogram()+
  jtools::theme_apa()+
  ggtitle("Duration by phase","units = log [seconds]")

Same plots but without “total surgery”.

Code 
dur %>% 
  pivot_longer(everything()) %>% 
  rename(phase = name, duration = value) %>% 
  filter(phase!="total surgery") %>% 
  ggplot(aes(x=duration))+
  facet_wrap(~phase,scales = "free",nrow = 2)+geom_histogram()+
  jtools::theme_apa()+
  ggtitle("Duration by phase","units = [seconds]")

Code 
dur %>% 
  pivot_longer(everything()) %>% 
  rename(phase = name, duration = value) %>% 
  filter(phase!="total surgery") %>% 
  ggplot(aes(x=log(duration)))+
  facet_wrap(~phase,scales = "free",nrow = 2)+geom_histogram()+
  jtools::theme_apa()+
    ggtitle("Duration by phase","units = log [seconds]")

Code 
pp <- 
dur %>% 
  pivot_longer(everything()) %>% 
  rename(phase = name, duration = value) %>% 
  filter(phase!="total surgery") %>% 
 ggplot(aes(x=log(duration)))+
facet_wrap(~phase,scales = "free")+geom_histogram()+
  ggtitle("Duration by phase","units = log [seconds]")

pp+  jtools::theme_apa()

By surgeon

We calculated summary statistics across operational phases for each surgeon.

  • A statistically significant differnce between median durations was found for every phase.
    • We use the Kruskal-Wallis rank-sum test (generalization of Mann-Whitney U-test) because duration data are non-normal.
    • The small p-values indicate rejection of \(H_0\) - median durations are equal across surgeons (for a given phase); We conclude that at least one surgeon differs from the others; To verify due to which surgeon this rejection occured, we conduct a post-hoc analysis.
Code 
dur %>% 
    add_column(surgeon = lean$surgeon) %>% 

  tbl_summary(by = surgeon,
                  type = all_continuous() ~ "continuous2",
    statistic = all_continuous() ~ c("{mean} ({sd})", "{median} ({p25}, {p75})", "{min}, {max}")) %>% 
  add_overall() %>% 
  add_p()
Characteristic Overall, N = 200 Surgeon1, N = 50 Surgeon2, N = 50 Surgeon3, N = 50 Surgeon4, N = 50 p-value1
surgical ports <0.001
    Mean (SD) 97 (44) 58 (23) 96 (35) 105 (39) 130 (45)
    Median (IQR) 91 (72, 111) 51 (43, 71) 88 (79, 99) 99 (81, 114) 116 (100, 143)
    Range 34, 271 34, 162 61, 238 57, 262 71, 271
capsulorhexis <0.001
    Mean (SD) 74 (32) 53 (21) 65 (14) 85 (43) 91 (26)
    Median (IQR) 67 (55, 84) 53 (42, 56) 62 (53, 73) 72 (66, 94) 85 (73, 103)
    Range 33, 336 33, 178 47, 110 50, 336 46, 175
hydrodissection and hydrodelineation <0.001
    Mean (SD) 71 (46) 35 (11) 57 (18) 91 (48) 101 (57)
    Median (IQR) 59 (42, 85) 32 (27, 42) 54 (47, 62) 78 (59, 115) 86 (72, 114)
    Range 23, 376 23, 77 32, 130 31, 286 31, 376
phacoemulsification <0.001
    Mean (SD) 282 (141) 207 (73) 222 (53) 265 (96) 434 (174)
    Median (IQR) 243 (191, 334) 192 (175, 216) 215 (179, 262) 242 (197, 319) 377 (319, 497)
    Range 87, 1,072 87, 520 143, 343 148, 688 225, 1,072
cortex aspiration <0.001
    Mean (SD) 140 (77) 96 (41) 114 (48) 137 (63) 213 (89)
    Median (IQR) 125 (84, 173) 82 (72, 101) 103 (78, 145) 124 (95, 154) 192 (154, 240)
    Range 43, 578 50, 225 49, 273 43, 322 117, 578
iol preparation <0.001
    Mean (SD) 92 (51) 57 (49) 84 (22) 106 (45) 121 (59)
    Median (IQR) 87 (64, 109) 46 (24, 69) 84 (69, 98) 99 (85, 109) 111 (88, 139)
    Range 13, 376 13, 260 34, 138 55, 358 37, 376
iol implantation <0.001
    Mean (SD) 70 (54) 59 (61) 43 (26) 59 (33) 117 (57)
    Median (IQR) 55 (40, 80) 50 (33, 65) 40 (27, 47) 54 (45, 58) 105 (86, 135)
    Range 15, 450 15, 450 19, 160 34, 206 34, 402
hydrating the ports and injection of antibiotic <0.001
    Mean (SD) 78 (70) 40 (23) 117 (70) 45 (30) 111 (90)
    Median (IQR) 57 (34, 93) 32 (22, 47) 98 (69, 139) 36 (29, 49) 84 (63, 126)
    Range 9, 521 9, 96 44, 393 14, 205 26, 521
total surgery <0.001
    Mean (SD) 905 (342) 606 (190) 799 (131) 895 (190) 1,320 (329)
    Median (IQR) 829 (659, 1,055) 557 (517, 635) 791 (699, 902) 873 (779, 979) 1,257 (1,109, 1,383)
    Range 432, 2,268 432, 1,627 588, 1,052 625, 1,545 778, 2,268
1 Kruskal-Wallis rank sum test

Almost all surgeons have unique means that differ from all the rest. Only for the top “ns” variables - you cannot say they differ.

Code 
dur %>% 
    add_column(surgeon = lean$surgeon) %>% 
  pivot_longer(cols = -surgeon,names_to = "phase",values_to =  "duration") %>% 
  group_by(phase) %>% 
  rstatix::pairwise_wilcox_test(duration ~ surgeon) %>% arrange(desc(p.adj)) %>%
  select(phase,group1,group2,p_value = p.adj.signif) %>% 
  pander::pander()
phase group1 group2 p_value
hydrating the ports and injection of antibiotic Surgeon1 Surgeon3 ns
hydrating the ports and injection of antibiotic Surgeon2 Surgeon4 ns
iol implantation Surgeon1 Surgeon3 ns
hydrodissection and hydrodelineation Surgeon3 Surgeon4 ns
surgical ports Surgeon2 Surgeon3 ns
iol preparation Surgeon3 Surgeon4 ns
phacoemulsification Surgeon1 Surgeon2 ns
cortex aspiration Surgeon1 Surgeon2 *
cortex aspiration Surgeon2 Surgeon3 *
phacoemulsification Surgeon2 Surgeon3 *
iol implantation Surgeon1 Surgeon2 *
capsulorhexis Surgeon3 Surgeon4 *
total surgery Surgeon2 Surgeon3 *
phacoemulsification Surgeon1 Surgeon3 **
surgical ports Surgeon3 Surgeon4 **
iol preparation Surgeon2 Surgeon3 ***
iol preparation Surgeon2 Surgeon4 ***
capsulorhexis Surgeon2 Surgeon3 ***
cortex aspiration Surgeon1 Surgeon3 ****
capsulorhexis Surgeon1 Surgeon2 ****
hydrodissection and hydrodelineation Surgeon2 Surgeon3 ****
iol preparation Surgeon1 Surgeon2 ****
iol implantation Surgeon2 Surgeon3 ****
surgical ports Surgeon2 Surgeon4 ****
cortex aspiration Surgeon3 Surgeon4 ****
capsulorhexis Surgeon2 Surgeon4 ****
phacoemulsification Surgeon3 Surgeon4 ****
hydrodissection and hydrodelineation Surgeon2 Surgeon4 ****
iol preparation Surgeon1 Surgeon4 ****
iol preparation Surgeon1 Surgeon3 ****
hydrating the ports and injection of antibiotic Surgeon3 Surgeon4 ****
hydrating the ports and injection of antibiotic Surgeon1 Surgeon4 ****
cortex aspiration Surgeon2 Surgeon4 ****
total surgery Surgeon1 Surgeon2 ****
surgical ports Surgeon1 Surgeon2 ****
hydrodissection and hydrodelineation Surgeon1 Surgeon2 ****
iol implantation Surgeon1 Surgeon4 ****
iol implantation Surgeon3 Surgeon4 ****
surgical ports Surgeon1 Surgeon3 ****
capsulorhexis Surgeon1 Surgeon3 ****
total surgery Surgeon3 Surgeon4 ****
hydrating the ports and injection of antibiotic Surgeon1 Surgeon2 ****
hydrating the ports and injection of antibiotic Surgeon2 Surgeon3 ****
cortex aspiration Surgeon1 Surgeon4 ****
total surgery Surgeon1 Surgeon3 ****
capsulorhexis Surgeon1 Surgeon4 ****
iol implantation Surgeon2 Surgeon4 ****
phacoemulsification Surgeon2 Surgeon4 ****
phacoemulsification Surgeon1 Surgeon4 ****
hydrodissection and hydrodelineation Surgeon1 Surgeon3 ****
surgical ports Surgeon1 Surgeon4 ****
hydrodissection and hydrodelineation Surgeon1 Surgeon4 ****
total surgery Surgeon2 Surgeon4 ****
total surgery Surgeon1 Surgeon4 ****

This fabulous plot summarizes your work graphically.

Code 
require(elucidate)
dur %>% 
  add_column(surgeon = lean$surgeon) %>% 
  pivot_longer(-surgeon) %>% 
  rename(phase = name, duration = value) -> ds

ds$phase <- fct(ds$phase)
ds$surgeon <- fct(ds$surgeon)
ds$phase <- fct_recode(ds$phase,
             "Surgical ports" = "surgical ports",
             "Capsulorhexis" = "capsulorhexis",
             "Hydrodissection & hydrodelineation" = "hydrodissection and hydrodelineation",
             "Phacoemulsification" = "phacoemulsification",
             "Cortex aspiration" = "cortex aspiration",
             "Cortex aspiration" = "iol preparation",
             "IOL implantation" = "iol preparation",
             "Hydrating the ports & injection of antibiotic" = "hydrating the ports and injection of antibiotic",
             "Total surgery" = "total surgery") 
ds$phase <- fct_inorder(ds$phase)
# ds$phase %>% fct_relevel(levels =c("Surgical ports",
#                            "Capsulorhexis",
#                            "Hydrodissection & hydrodelineation",
#                            "Phacoemulsification",
#                           "Cortex aspiration",
#                           "Cortex aspiration",
#                            "IOL implantation",
#                           "Hydrating the ports & injection of antibiotic",
#                           "Total surgery"))
elucidate::plot_raincloud(data = ds %>% filter(phase != "total surgery"),x = surgeon,y=duration,coord_flip = T,fill_var = surgeon,
                          fill_var_order = rev(levels(ds$surgeon)),xlab = "",title = "Operation phase durations [sec.]",omit_legend = F,legend_position = "left")+#,facet_var = "phase",
                          #facet_var_order =levels(ds$phase))+
   
  facet_wrap(~phase,nrow = 2,scales = "free_x")#+

Code 
  theme(axis.text.y = element_blank()) +
  theme(panel.border = element_blank(), panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) +
  theme(legend.title = element_blank())
List of 6
 $ axis.text.y     : list()
  ..- attr(*, "class")= chr [1:2] "element_blank" "element"
 $ axis.line       :List of 6
  ..$ colour       : chr "black"
  ..$ linewidth    : NULL
  ..$ linetype     : NULL
  ..$ lineend      : NULL
  ..$ arrow        : logi FALSE
  ..$ inherit.blank: logi FALSE
  ..- attr(*, "class")= chr [1:2] "element_line" "element"
 $ panel.border    : list()
  ..- attr(*, "class")= chr [1:2] "element_blank" "element"
 $ panel.grid.major: list()
  ..- attr(*, "class")= chr [1:2] "element_blank" "element"
 $ panel.grid.minor: list()
  ..- attr(*, "class")= chr [1:2] "element_blank" "element"
 $ legend.title    : list()
  ..- attr(*, "class")= chr [1:2] "element_blank" "element"
 - attr(*, "class")= chr [1:2] "theme" "gg"
 - attr(*, "complete")= logi FALSE
 - attr(*, "validate")= logi TRUE
Code 
dur %>% 
  add_column(surgeon = lean$surgeon) %>% 
  pivot_longer(-surgeon) %>% 
  rename(phase = name, duration = value) %>% 
  filter(phase == "total surgery") -> td
td$surgeon <- fct_rev(td$surgeon)
elucidate::plot_raincloud(
  data = td,
  y = duration,
  x = surgeon,x_var_order = levels(td$surgeon),
  box_plot = T,
  point_size = 2,
  coord_flip = T,
  box_side = "l",
  box_colour = "green",
  box_nudge = .3,
  box_outlier_shape = 14,
  box_outlier_size = 0,
  box_line_size = .5,
  title = "Overall surgery duration by surgeon")+
  ylab("Overall duration [sec.]") +
  xlab('') +
  theme(
    panel.border = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.line = element_line(colour = "black")
  )

Phase duration by surgeon :

Code 
surgeon_number <- 1
  surg_level <- paste0("Surgeon",surgeon_number,collapse = "")
  tit <- paste0("Durations for surgeon ", surgeon_number, collapse = "")
  ds %>%
  filter(phase!= "Total surgery") %>% 
  filter(surgeon==surg_level) %>% 
  
  elucidate::plot_raincloud(x = phase, y = duration,coord_flip = T,box_plot = T,
  point_size = 2,
  box_side = "l",
  box_colour = "green",
  box_nudge = .3,box_line_size = 0.5,
  box_outlier_fill = NA,box_outlier_size = 0,
  title =tit)

Code 
  surgeon_number <- 2
  surg_level <- paste0("Surgeon",surgeon_number,collapse = "")
  tit <- paste0("Durations for surgeon ", surgeon_number, collapse = "")
  ds %>%
  filter(phase!= "Total surgery") %>% 
  filter(surgeon==surg_level) %>% 
  
  elucidate::plot_raincloud(x = phase, y = duration,coord_flip = T,box_plot = T,
  point_size = 2,
  box_side = "l",
  box_colour = "green",
  box_nudge = .3,box_line_size = 0.5,
  box_outlier_fill = NA,box_outlier_size = 0,
  title =tit)

Code 
  surgeon_number <- 3
  surg_level <- paste0("Surgeon",surgeon_number,collapse = "")
  tit <- paste0("Durations for surgeon ", surgeon_number, collapse = "")
  ds %>%
  filter(phase!= "Total surgery") %>% 
  filter(surgeon==surg_level) %>% 
  
  elucidate::plot_raincloud(x = phase, y = duration,coord_flip = T,box_plot = T,
  point_size = 2,
  box_side = "l",
  box_colour = "green",
  box_nudge = .3,box_line_size = 0.5,
  box_outlier_fill = NA,box_outlier_size = 0,
  title =tit)

Code 
  surgeon_number <- 4
  surg_level <- paste0("Surgeon",surgeon_number,collapse = "")
  tit <- paste0("Durations for surgeon ", surgeon_number, collapse = "")
  ds %>%
  filter(phase!= "Total surgery") %>% 
  filter(surgeon==surg_level) %>% 
  
  elucidate::plot_raincloud(x = phase, y = duration,coord_flip = T,box_plot = T,
  point_size = 2,
  box_side = "l",
  box_colour = "green",
  box_nudge = .3,box_line_size = 0.5,
  box_outlier_fill = NA,box_outlier_size = 0,
  title =tit)

By complexity

  • Complexity data was available only for RMC surgeons (3 and 4), thus the analyses below have limited generelization power.
Code 
dur %>% 
  add_column(surgeon = lean$surgeon,
             complexity = lean$complexity) %>% 
  mutate(complexity = fct_explicit_na(factor(complexity))) %>% 
  pivot_longer(-c(surgeon,complexity)) %>% 
  rename(phase = name, duration = value) %>% 
  filter(phase != "total surgery") -> tc

lean %>% 
  tbl_cross(row = surgeon, col = complexity) %>% 
  modify_caption("Surgery complexity * Surgeon crosstabulation")
Surgery complexity * Surgeon crosstabulation
complexity Total
1 2 3 Unknown
surgeon
    Surgeon1 0 0 0 50 50
    Surgeon2 0 0 0 50 50
    Surgeon3 20 21 4 5 50
    Surgeon4 14 22 8 6 50
Total 34 43 12 111 200

  • The following tables reveal which phases’ durations are associated with the complexity.

  • Analyses with and without the missing observations (complexity unknown in 111 of 198) lead to different conclusions.

  • OMITTING the missing observations (first table) we have:

  • total duration (p=.015) and phacoemulsification (p=.017) are longer for complexity = 3.

  • surgical ports (p=.057) and iol_implantation (p = .078) are marginally significant.

  • The remaining phases’ durations are not associated with complexity.

Code 
dur %>% 
  add_column(
             complexity = lean$complexity) %>% 
  tbl_summary(by = complexity,
                  type = all_continuous() ~ "continuous2",
    statistic = all_continuous() ~ c("{mean} ({sd})", "{median} ({p25}, {p75})", "{min}, {max}")) %>% 
  add_overall() %>% 
  add_p() %>% 
  sort_p()
Characteristic Overall, N = 89 1, N = 34 2, N = 43 3, N = 12 p-value1
total surgery 0.004
    Mean (SD) 1,102 (333) 1,008 (257) 1,107 (360) 1,351 (315)
    Median (IQR) 1,031 (875, 1,257) 963 (830, 1,202) 1,021 (882, 1,176) 1,272 (1,223, 1,479)
    Range 625, 2,268 625, 1,620 709, 2,268 830, 2,164
phacoemulsification 0.011
    Mean (SD) 345 (160) 302 (126) 350 (164) 450 (193)
    Median (IQR) 319 (242, 377) 273 (200, 372) 314 (260, 364) 365 (345, 506)
    Range 148, 1,072 162, 688 148, 1,072 294, 951
surgical ports 0.014
    Mean (SD) 117 (45) 104 (34) 117 (37) 157 (73)
    Median (IQR) 106 (90, 125) 99 (84, 113) 109 (94, 125) 134 (109, 209)
    Range 57, 271 57, 229 62, 240 70, 271
iol implantation 0.037
    Mean (SD) 86 (47) 74 (36) 87 (48) 121 (54)
    Median (IQR) 71 (52, 110) 59 (46, 96) 61 (55, 108) 130 (87, 157)
    Range 34, 208 34, 180 36, 208 34, 206
cortex aspiration 0.2
    Mean (SD) 175 (86) 170 (67) 173 (105) 196 (60)
    Median (IQR) 155 (118, 196) 158 (127, 193) 150 (114, 192) 203 (150, 239)
    Range 65, 578 81, 322 65, 578 95, 281
capsulorhexis 0.2
    Mean (SD) 88 (37) 95 (50) 81 (25) 91 (29)
    Median (IQR) 81 (67, 101) 83 (71, 104) 73 (66, 93) 92 (71, 102)
    Range 46, 336 57, 336 46, 175 50, 163
iol preparation 0.2
    Mean (SD) 112 (47) 106 (33) 108 (38) 147 (87)
    Median (IQR) 103 (85, 129) 100 (85, 118) 102 (85, 124) 125 (94, 166)
    Range 37, 358 50, 200 37, 222 55, 358
hydrodissection and hydrodelineation 0.6
    Mean (SD) 98 (55) 88 (40) 104 (66) 102 (46)
    Median (IQR) 85 (63, 114) 84 (55, 111) 84 (66, 118) 95 (65, 120)
    Range 31, 376 31, 173 31, 376 54, 213
hydrating the ports and injection of antibiotic >0.9
    Mean (SD) 78 (78) 69 (45) 87 (99) 71 (68)
    Median (IQR) 51 (34, 88) 63 (37, 89) 51 (34, 93) 50 (34, 74)
    Range 14, 521 14, 200 18, 521 23, 276
1 Kruskal-Wallis rank sum test
  • Keeping the missing observations and treating them as a separate class:
    • All phases appear to have significant associations of complexity and duration.
    • Note how the median and mean of the missing group are smaller than ALL other three complexity levels. It is unclear whether this stems from the compte
Code 
dur %>% 
  add_column(
             complexity = lean$complexity) %>% 
  mutate(complexity = fct_explicit_na(factor(complexity))) %>% 
  tbl_summary(by = complexity,
                  type = all_continuous() ~ "continuous2",
    statistic = all_continuous() ~ c("{mean} ({sd})", "{median} ({p25}, {p75})", "{min}, {max}")) %>% 
  add_overall() %>% 
  add_p() %>% 
  sort_p()
Characteristic Overall, N = 200 1, N = 34 2, N = 43 3, N = 12 (Missing), N = 111 p-value1
total surgery <0.001
    Mean (SD) 905 (342) 1,008 (257) 1,107 (360) 1,351 (315) 747 (259)
    Median (IQR) 829 (659, 1,055) 963 (830, 1,202) 1,021 (882, 1,176) 1,272 (1,223, 1,479) 697 (580, 843)
    Range 432, 2,268 625, 1,620 709, 2,268 830, 2,164 432, 2,136
hydrodissection and hydrodelineation <0.001
    Mean (SD) 71 (46) 88 (40) 104 (66) 102 (46) 50 (22)
    Median (IQR) 59 (42, 85) 84 (55, 111) 84 (66, 118) 95 (65, 120) 47 (32, 60)
    Range 23, 376 31, 173 31, 376 54, 213 23, 130
capsulorhexis <0.001
    Mean (SD) 74 (32) 95 (50) 81 (25) 91 (29) 62 (21)
    Median (IQR) 67 (55, 84) 83 (71, 104) 73 (66, 93) 92 (71, 102) 56 (50, 69)
    Range 33, 336 57, 336 46, 175 50, 163 33, 178
phacoemulsification <0.001
    Mean (SD) 282 (141) 302 (126) 350 (164) 450 (193) 232 (99)
    Median (IQR) 243 (191, 334) 273 (200, 372) 314 (260, 364) 365 (345, 506) 205 (178, 257)
    Range 87, 1,072 162, 688 148, 1,072 294, 951 87, 764
surgical ports <0.001
    Mean (SD) 97 (44) 104 (34) 117 (37) 157 (73) 82 (37)
    Median (IQR) 91 (72, 111) 99 (84, 113) 109 (94, 125) 134 (109, 209) 77 (54, 95)
    Range 34, 271 57, 229 62, 240 70, 271 34, 238
cortex aspiration <0.001
    Mean (SD) 140 (77) 170 (67) 173 (105) 196 (60) 112 (54)
    Median (IQR) 125 (84, 173) 158 (127, 193) 150 (114, 192) 203 (150, 239) 92 (77, 138)
    Range 43, 578 81, 322 65, 578 95, 281 43, 290
iol implantation <0.001
    Mean (SD) 70 (54) 74 (36) 87 (48) 121 (54) 56 (56)
    Median (IQR) 55 (40, 80) 59 (46, 96) 61 (55, 108) 130 (87, 157) 45 (33, 61)
    Range 15, 450 34, 180 36, 208 34, 206 15, 450
iol preparation <0.001
    Mean (SD) 92 (51) 106 (33) 108 (38) 147 (87) 76 (49)
    Median (IQR) 87 (64, 109) 100 (85, 118) 102 (85, 124) 125 (94, 166) 72 (50, 94)
    Range 13, 376 50, 200 37, 222 55, 358 13, 376
hydrating the ports and injection of antibiotic >0.9
    Mean (SD) 78 (70) 69 (45) 87 (99) 71 (68) 78 (63)
    Median (IQR) 57 (34, 93) 63 (37, 89) 51 (34, 93) 50 (34, 74) 62 (34, 100)
    Range 9, 521 14, 200 18, 521 23, 276 9, 393
1 Kruskal-Wallis rank sum test
Code 
td <- 
  td %>% add_column(complexity = lean$complexity)

elucidate::plot_raincloud(
  data = td,
  y = duration,
  x = complexity,
  box_plot = T,
  point_size = 2,
  coord_flip = T,
  box_side = "l",
  box_colour = "green",box_outlier_size = 0,
  box_nudge = .3,
  box_line_size = .5,
  title = "Overall surgery duration by complexity"
) + ylab("Overall duration [sec.]") +
  xlab('') +
  theme(
    panel.border = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.line = element_line(colour = "black")
  )

Code 
elucidate::plot_raincloud(
  data = dur %>% 
  add_column(surgeon = lean$surgeon, 
             complexity = lean$complexity,
             duration = lean$total_surgery_duration),
             y = duration,
  x = complexity,
  facet_var = surgeon,
  box_plot = T,
  point_size = 2,
  coord_flip = T,
  box_side = "l",
  box_colour = "green",box_outlier_size = 0,
  box_nudge = .3,
  box_line_size = .5,
  title = "Overall surgery duration by complexity & Surgeon"
) + ylab("Overall duration [sec.]") +
  xlab('') +
  theme(
    panel.border = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.line = element_line(colour = "black")
  )

By medical center

Code 
tbl_cross(lean,row = or_location,col = surgeon)
surgeon Total
Surgeon1 Surgeon2 Surgeon3 Surgeon4
or_location
    Maayan 50 50 0 0 100
    Rambam 0 0 50 50 100
Total 50 50 50 50 200

You know why RMC is slower than Maayan right?

Maybe it’s the because of surgeon 4 who works here ?

Code 
dur %>% 
  add_column(hospital = lean$or_location) %>% 
  tbl_summary(by = hospital,
              type = all_continuous() ~ "continuous2",
    statistic = all_continuous() ~ c("{mean} ({sd})", "{median} ({p25}, {p75})", "{min}, {max}")) %>% 
  add_p()
Characteristic Maayan, N = 100 Rambam, N = 100 p-value1
surgical ports <0.001
    Mean (SD) 77 (35) 118 (44)
    Median (IQR) 75 (51, 91) 109 (91, 127)
    Range 34, 238 57, 271
capsulorhexis <0.001
    Mean (SD) 59 (19) 88 (36)
    Median (IQR) 55 (49, 65) 81 (68, 101)
    Range 33, 178 46, 336
hydrodissection and hydrodelineation <0.001
    Mean (SD) 46 (19) 96 (52)
    Median (IQR) 45 (32, 58) 84 (63, 114)
    Range 23, 130 31, 376
phacoemulsification <0.001
    Mean (SD) 215 (64) 350 (164)
    Median (IQR) 203 (176, 245) 319 (226, 396)
    Range 87, 520 148, 1,072
cortex aspiration <0.001
    Mean (SD) 105 (45) 175 (85)
    Median (IQR) 87 (74, 128) 154 (118, 214)
    Range 49, 273 43, 578
iol preparation <0.001
    Mean (SD) 71 (40) 114 (53)
    Median (IQR) 69 (42, 91) 103 (85, 127)
    Range 13, 260 37, 376
iol implantation <0.001
    Mean (SD) 51 (47) 88 (55)
    Median (IQR) 43 (31, 57) 69 (53, 108)
    Range 15, 450 34, 402
hydrating the ports and injection of antibiotic 0.7
    Mean (SD) 79 (64) 78 (75)
    Median (IQR) 61 (32, 98) 54 (34, 89)
    Range 9, 393 14, 521
total surgery <0.001
    Mean (SD) 703 (189) 1,108 (342)
    Median (IQR) 683 (558, 810) 1,033 (869, 1,261)
    Range 432, 1,627 625, 2,268
1 Wilcoxon rank sum test

Even if we remove surgeon 4 and use only surgeon 3, the difference is still in favor of maayan MC as seen in the table below. However, surgeon 3 was the second slowest on average.

To summarize - the between-surgeon and between-medical center effects cannot be distinguished because the surgeon is nested in a single medical center.

Code 
dur %>% 
   add_column(surgeon = lean$surgeon) %>% 
  add_column(hospital = lean$or_location) %>% 
    filter(surgeon != "Surgeon4") %>% 
  tbl_summary(by = hospital,
              type = all_continuous() ~ "continuous2",
    statistic = all_continuous() ~ c("{mean} ({sd})", "{median} ({p25}, {p75})", "{min}, {max}")) %>% 
  add_p()
Characteristic Maayan, N = 1001 Rambam, N = 501 p-value2
surgical ports <0.001
    Mean (SD) 77 (35) 105 (39)
    Median (IQR) 75 (51, 91) 99 (81, 114)
    Range 34, 238 57, 262
capsulorhexis <0.001
    Mean (SD) 59 (19) 85 (43)
    Median (IQR) 55 (49, 65) 72 (66, 94)
    Range 33, 178 50, 336
hydrodissection and hydrodelineation <0.001
    Mean (SD) 46 (19) 91 (48)
    Median (IQR) 45 (32, 58) 78 (59, 115)
    Range 23, 130 31, 286
phacoemulsification <0.001
    Mean (SD) 215 (64) 265 (96)
    Median (IQR) 203 (176, 245) 242 (197, 319)
    Range 87, 520 148, 688
cortex aspiration <0.001
    Mean (SD) 105 (45) 137 (63)
    Median (IQR) 87 (74, 128) 124 (95, 154)
    Range 49, 273 43, 322
iol preparation <0.001
    Mean (SD) 71 (40) 106 (45)
    Median (IQR) 69 (42, 91) 99 (85, 109)
    Range 13, 260 55, 358
iol implantation <0.001
    Mean (SD) 51 (47) 59 (33)
    Median (IQR) 43 (31, 57) 54 (45, 58)
    Range 15, 450 34, 206
hydrating the ports and injection of antibiotic <0.001
    Mean (SD) 79 (64) 45 (30)
    Median (IQR) 61 (32, 98) 36 (29, 49)
    Range 9, 393 14, 205
total surgery <0.001
    Mean (SD) 703 (189) 895 (190)
    Median (IQR) 683 (558, 810) 873 (779, 979)
    Range 432, 1,627 625, 1,545
surgeon <0.001
    Surgeon1 50 (50%) 0 (0%)
    Surgeon2 50 (50%) 0 (0%)
    Surgeon3 0 (0%) 50 (100%)
1 n (%)
2 Wilcoxon rank sum test; Pearson’s Chi-squared test

NEW - Personal surgeon development

I propose two approaches to quantify a surgeon’s performance.

  1. “Consistency” using CV - measure the coefficient of variation. Smaller values indicate a more consistent performance, with similar durations.

  2. “Efficiency” using frequency - how many surgeries/procedure can the surgeon cram into one time unit? The more, the better. I define frequency as: \(f = 3600/Y\)

where \(Y\) is the duration of a phase. The units are how-many-per-hour would the surgeon be able to do if he worked in the current rate.

Approach 1

Code 
dat <- read_csv("LEAN_final_210922.csv") 

cv <- function(x,...) return(sd(x) / mean(x))
CV <- Vectorize(cv)
dat %>%
  select(-zepto, -surgery_date, -file_name) %>%
  tbl_summary(
    include = ends_with("duration"),
    statistic = list(all_continuous() ~ "{mean}± {sd} || {cv}"),
    digits = list(all_continuous() ~ c(1,2,2))
  )
Characteristic N = 2001
surgical_ports_duration 97.4± 44.41 || 0.46
capsulorhexis_duration 73.6± 31.82 || 0.43
hydrodissection_and_hydrodelineation_duration 71.2± 46.48 || 0.65
phacoemulsification_duration 282.1± 141.18 || 0.50
cortex_aspiration_duration 140.0± 76.74 || 0.55
iol_preparation_duration 92.2± 51.48 || 0.56
iol_implantation_duration 69.6± 54.21 || 0.78
hydrating_the_ports_and_injection_of_antibiotic_duration 78.2± 69.58 || 0.89
total_surgery_duration 905.2± 342.43 || 0.38
1 Mean± SD || cv
Code 
dat %>%
  select(-zepto, -surgery_date, -file_name) %>%
  tbl_summary(
    by = surgeon,
    include = ends_with("duration"),
    statistic = list(all_continuous() ~ "{mean}± {sd} || {cv}"),
    digits = list(all_continuous() ~ c(1,2,2))
  ) %>% 
  add_overall()
Characteristic Overall, N = 2001 1, N = 501 2, N = 501 3, N = 501 4, N = 501
surgical_ports_duration 97.4± 44.41 || 0.46 58.4± 22.91 || 0.39 96.3± 34.62 || 0.36 105.4± 39.50 || 0.37 129.6± 45.19 || 0.35
capsulorhexis_duration 73.6± 31.82 || 0.43 53.5± 21.12 || 0.40 64.8± 14.37 || 0.22 85.4± 43.34 || 0.51 90.5± 25.59 || 0.28
hydrodissection_and_hydrodelineation_duration 71.2± 46.48 || 0.65 35.3± 11.41 || 0.32 57.3± 17.94 || 0.31 91.4± 48.15 || 0.53 100.6± 56.50 || 0.56
phacoemulsification_duration 282.1± 141.18 || 0.50 207.1± 73.13 || 0.35 222.2± 52.91 || 0.24 265.0± 96.12 || 0.36 434.1± 174.04 || 0.40
cortex_aspiration_duration 140.0± 76.74 || 0.55 95.5± 40.58 || 0.42 114.2± 48.13 || 0.42 137.4± 62.90 || 0.46 212.8± 88.76 || 0.42
iol_preparation_duration 92.2± 51.48 || 0.56 57.2± 49.05 || 0.86 84.4± 22.08 || 0.26 105.9± 45.30 || 0.43 121.3± 58.76 || 0.48
iol_implantation_duration 69.6± 54.21 || 0.78 58.9± 60.70 || 1.03 43.3± 26.08 || 0.60 59.4± 32.91 || 0.55 116.8± 57.36 || 0.49
hydrating_the_ports_and_injection_of_antibiotic_duration 78.2± 69.58 || 0.89 40.3± 23.02 || 0.57 116.7± 69.61 || 0.60 44.8± 30.38 || 0.68 111.1± 90.27 || 0.81
total_surgery_duration 905.2± 342.43 || 0.38 606.2± 189.89 || 0.31 799.2± 131.04 || 0.16 894.8± 189.51 || 0.21 1,320.4± 329.10 || 0.25
1 Mean± SD || cv

Table of CV by surgeon * duration with 95% percentile Bootstrap confidence intervals with \(B=10,000\) bootstrap resamples.

Code 
sdur <- 
dur %>% 
  add_column(surgeon = dat$surgeon) %>% 
  pivot_longer(-surgeon,names_to = "phase",values_to = "duration")

cv0 <- 
sdur %>% 
  group_by(surgeon,phase) %>% 
  summarize(cv = sd(duration)/mean(duration)) 

cvboot <- replicate(
  1e3,
  sdur %>%
    
    group_by(surgeon, phase) %>%
    sample_frac(replace = TRUE) %>%
    summarize(cv = sd(duration)/mean(duration)) %>%
    ungroup() %>%
    pull(cv)
  
)

cvboot %>% 
  t() %>% 
  apply(., 2, quantile,c(.025,0.975)) %>% 
  t() %>% 
  as_tibble() %>% 
  set_names(c("lower","upper")) %>% 
  add_column(cv0) %>%
  mutate(across(all_of(c("lower","upper","cv")),~round(.x,2))) %>% 
  mutate(
    CV = str_glue("{cv} [{lower}-{upper}]")
  ) %>% 
  select(surgeon,phase,CV) %>% 
  pivot_wider(id_cols = c("surgeon"),values_from = ("CV"),names_from =  c("phase")) -> cv_tab


cv_tab <- 
sjmisc::rotate_df(cv_tab,cn = T) %>% 
    set_names(paste0("Surgeon ",1:4))


pander::pander(cv_tab)
Table continues below
  Surgeon 1 Surgeon 2
capsulorhexis 0.4 [0.18-0.58] 0.22 [0.17-0.26]
cortex aspiration 0.42 [0.33-0.49] 0.42 [0.33-0.5]
hydrating the ports and injection of antibiotic 0.57 [0.47-0.64] 0.6 [0.43-0.73]
hydrodissection and hydrodelineation 0.32 [0.24-0.39] 0.31 [0.2-0.4]
iol implantation 1.03 [0.37-1.37] 0.6 [0.29-0.75]
iol preparation 0.86 [0.57-1.04] 0.26 [0.21-0.31]
phacoemulsification 0.35 [0.21-0.45] 0.24 [0.2-0.27]
surgical ports 0.39 [0.27-0.51] 0.36 [0.17-0.46]
total surgery 0.31 [0.14-0.44] 0.16 [0.14-0.18]
  Surgeon 3 Surgeon 4
capsulorhexis 0.51 [0.22-0.7] 0.28 [0.21-0.34]
cortex aspiration 0.46 [0.35-0.53] 0.42 [0.28-0.51]
hydrating the ports and injection of antibiotic 0.68 [0.36-0.87] 0.81 [0.5-0.98]
hydrodissection and hydrodelineation 0.53 [0.4-0.64] 0.56 [0.34-0.72]
iol implantation 0.55 [0.2-0.69] 0.49 [0.3-0.64]
iol preparation 0.43 [0.18-0.6] 0.48 [0.34-0.61]
phacoemulsification 0.36 [0.26-0.46] 0.4 [0.28-0.48]
surgical ports 0.37 [0.24-0.47] 0.35 [0.26-0.41]
total surgery 0.21 [0.15-0.26] 0.25 [0.19-0.29]

Approach 2

The units of the following table are “procedures per hour”.

Code 
# dur %>% 
#   map_df(~3600/.x) -> dur_inv
dur_inv <- dur
dur_inv <- janitor::clean_names(dur_inv)


wsdi <- dur_inv %>% 
  add_column(surgeon = lean$surgeon)

wsdi %>% 
  tbl_summary(by = surgeon)
Characteristic Surgeon1, N = 501 Surgeon2, N = 501 Surgeon3, N = 501 Surgeon4, N = 501
surgical_ports 51 (43, 71) 88 (79, 99) 99 (81, 114) 116 (100, 143)
capsulorhexis 53 (42, 56) 62 (53, 73) 72 (66, 94) 85 (73, 103)
hydrodissection_and_hydrodelineation 32 (27, 42) 54 (47, 62) 78 (59, 115) 86 (72, 114)
phacoemulsification 192 (175, 216) 215 (179, 262) 242 (197, 319) 377 (319, 497)
cortex_aspiration 82 (72, 101) 103 (78, 145) 124 (95, 154) 192 (154, 240)
iol_preparation 46 (24, 69) 84 (69, 98) 99 (85, 109) 111 (88, 139)
iol_implantation 50 (33, 65) 40 (27, 47) 54 (45, 58) 105 (86, 135)
hydrating_the_ports_and_injection_of_antibiotic 32 (22, 47) 98 (69, 139) 36 (29, 49) 84 (63, 126)
total_surgery 557 (517, 635) 791 (699, 902) 873 (779, 979) 1,257 (1,109, 1,383)
1 Median (IQR)
Code 
elucidate::plot_var_all(wsdi,group_var = surgeon)

Footnotes

  1. The analyst is not sure how much freedom of choice surgeons have in this procedure.↩︎