SNAP Audit

1 Setup

1.1 Load in Data

snap <- fread("snap_audit.csv",
              na.strings = c(
                "Not available frominformation available",
                "Not available from information available",
                "Not clear from information available",
                "Unclear",
                "Not Clear from the information available",
                "",
                " ",
                "Not clear fro the information available",
                "Not clear from the information available",
                "NA",
                "na",
                "N/A"
                )) %>% 
  as_tibble() %>%
  janitor::clean_names() %>% 
  subset(select = -c(ip,
                     who_are_you)) %>% 
  mutate(
    submission_id = 1:n()
  ) 

## Warning in fread("snap_audit.csv", na.strings = c("Not available
## frominformation available", : na.strings[7]==" " consists only of whitespace,
## ignoring. strip.white==TRUE (default) and "" is present in na.strings, so any
## number of spaces in string columns will already be read as <NA>.

2 Explore outcomes - initial look

2.1 Plot missing

outcomes <- snap %>%
  subset(
    select = c(
      submission_id,
      at_follow_up_did_the_patient_describe_ongoing_pain,
      did_the_patient_describe_new_or_worsening_erectile_dysfunction_at_follow_up,
      at_follow_up_did_the_patient_describe_penile_curvature,
      at_follow_up_was_there_evidence_of_penile_wasting_or_shortening,
      were_there_any_early_complications_clavien_dindo_3
    )
  ) %>% 
  mutate(good_outcome = case_when(
    at_follow_up_did_the_patient_describe_ongoing_pain == "No" & 
      did_the_patient_describe_new_or_worsening_erectile_dysfunction_at_follow_up == "No" & 
      at_follow_up_did_the_patient_describe_penile_curvature == "No" & 
      at_follow_up_was_there_evidence_of_penile_wasting_or_shortening == "No" & 
      were_there_any_early_complications_clavien_dindo_3 == "No" ~ "Yes",
    at_follow_up_did_the_patient_describe_ongoing_pain == "Yes" |
      did_the_patient_describe_new_or_worsening_erectile_dysfunction_at_follow_up == "Yes" | 
      at_follow_up_did_the_patient_describe_penile_curvature == "Yes" |
      at_follow_up_was_there_evidence_of_penile_wasting_or_shortening == "Yes" |
      were_there_any_early_complications_clavien_dindo_3 == "Yes" ~ "No",
    TRUE ~ NA_character_
  ))

plot_missing(outcomes)

2.2 Pain

outcomes %>% 
  drop_na(at_follow_up_did_the_patient_describe_ongoing_pain) %>%
  filter(at_follow_up_did_the_patient_describe_ongoing_pain != "No follow up" &
           at_follow_up_did_the_patient_describe_ongoing_pain != "No Follow up" &
           at_follow_up_did_the_patient_describe_ongoing_pain != "Not from the information available") %>%
  tabyl(at_follow_up_did_the_patient_describe_ongoing_pain) %>%
  mutate(
    Pain = at_follow_up_did_the_patient_describe_ongoing_pain,
    Percent = round(percent * 100),
         .keep = "unused") %>%
  select(c(Pain,
             n,
             Percent)) %>%
  gt() %>%
  gt_theme_espn()

Pain	n	Percent
No	130	82
Yes	28	18

2.3 New/Worsening ED

outcomes %>% 
  drop_na(did_the_patient_describe_new_or_worsening_erectile_dysfunction_at_follow_up) %>%
  subset(did_the_patient_describe_new_or_worsening_erectile_dysfunction_at_follow_up != "No follow up") %>%
  tabyl(did_the_patient_describe_new_or_worsening_erectile_dysfunction_at_follow_up) %>%
  mutate(
    "New/Worsening ED" = did_the_patient_describe_new_or_worsening_erectile_dysfunction_at_follow_up,
    Percent = round(percent * 100),
         .keep = "unused") %>%
  select(c("New/Worsening ED",
             n,
             Percent)) %>%
  gt() %>%
  gt_theme_espn()

New/Worsening ED	n	Percent
No	113	74
Yes - Manageable with medication	23	15
Yes - Still bothersome despite medication	5	3
Yes - Unable to quantify from information available	11	7

2.4 Penile Curvature

outcomes %>% 
  drop_na(at_follow_up_did_the_patient_describe_penile_curvature) %>%
  subset(at_follow_up_did_the_patient_describe_penile_curvature != "No follow up") %>%
  tabyl(at_follow_up_did_the_patient_describe_penile_curvature) %>%
  mutate(
    "Penile Curvature" = at_follow_up_did_the_patient_describe_penile_curvature,
    Percent = round(percent * 100),
         .keep = "unused") %>%
  select(c("Penile Curvature",
             n,
             Percent)) %>%
  gt() %>%
  gt_theme_espn()

Penile Curvature	n	Percent
No	104	72
Yes - Able to have penetrative intercourse	35	24
Yes - Not clear if limiting penetrative intercourse	5	3

2.5 Penile Waisting / Shortening

outcomes %>% 
  drop_na(at_follow_up_was_there_evidence_of_penile_wasting_or_shortening) %>%
  subset(at_follow_up_was_there_evidence_of_penile_wasting_or_shortening != "No follow up") %>%
  tabyl(at_follow_up_was_there_evidence_of_penile_wasting_or_shortening) %>%
  mutate(
    "Penile Waisting / Shortening" = at_follow_up_was_there_evidence_of_penile_wasting_or_shortening,
    Percent = round(percent * 100),
         .keep = "unused") %>%
  select(c("Penile Waisting / Shortening",
             n,
             Percent)) %>%
  gt() %>%
  gt_theme_espn()

Penile Waisting / Shortening	n	Percent
No	114	87
Yes	17	13

2.6 Clavien Dindo >/=III

outcomes %>% 
  drop_na(were_there_any_early_complications_clavien_dindo_3) %>%
  tabyl(were_there_any_early_complications_clavien_dindo_3) %>%
  mutate(
    "Clavien Dindo >/= III" = were_there_any_early_complications_clavien_dindo_3,
    Percent = round(percent * 100),
         .keep = "unused") %>%
  select(c("Clavien Dindo >/= III",
             n,
             Percent)) %>%
  gt() %>%
  gt_theme_espn()

Clavien Dindo >/= III	n	Percent
No	245	98
Yes	5	2

2.7 Composite Good Outcome

outcomes %>% 
  drop_na(good_outcome) %>%
  tabyl(good_outcome) %>%
  mutate(
    "Composite Good Outcome" = good_outcome,
    Percent = round(percent * 100),
         .keep = "unused") %>%
  select(c("Composite Good Outcome",
             n,
             Percent)) %>%
  gt() %>%
  gt_theme_espn()

Composite Good Outcome	n	Percent
No	46	46
Yes	53	54

2.8 Define data

There are some dates that appear to be prior to audit collection?

snap <- snap %>%
  rowwise() %>%
  mutate(
    what_was_the_estimated_date_and_time_of_injury = case_when(
      str_detect(what_was_the_estimated_date_and_time_of_injury, ":") ~
        dmy_hm(what_was_the_estimated_date_and_time_of_injury),
      TRUE ~ dmy(what_was_the_estimated_date_and_time_of_injury) + hours(12)
    ),
    
    what_age_was_the_patient_at_the_time_of_the_injury = as.integer(what_age_was_the_patient_at_the_time_of_the_injury),
    mechanism_of_injury = as.factor(mechanism_of_injury),
    if_the_injury_was_sustained_during_sexual_intercourse_what_sexual_position_was_involved = as.factor(
      if_the_injury_was_sustained_during_sexual_intercourse_what_sexual_position_was_involved
    ),
    did_the_patient_have_a_history_of_a_prior_penile_fracture = as.factor(did_the_patient_have_a_history_of_a_prior_penile_fracture),
    when_was_the_patient_first_seen_in_hospital_by_a_healthcare_professional = case_when(
      str_detect(
        when_was_the_patient_first_seen_in_hospital_by_a_healthcare_professional,
        ":"
      ) ~
        dmy_hm(
          when_was_the_patient_first_seen_in_hospital_by_a_healthcare_professional
        ),
      TRUE ~ dmy(
        when_was_the_patient_first_seen_in_hospital_by_a_healthcare_professional
      ) + hours(12)
    ),
    was_pain_present_at_presentation = as.factor(was_pain_present_at_presentation),
    was_there_immediate_detumescence_loss_of_erection = as.factor(was_there_immediate_detumescence_loss_of_erection),
    was_an_audible_pop_reported_by_the_patient = as.factor(was_an_audible_pop_reported_by_the_patient),
    was_penile_bruising_present = as.factor(was_penile_bruising_present),
    did_bruising_extend_over_the_scrotum_lower_abdomen_or_perineum_answer_more_than_one_if_necessary = as.factor(
      did_bruising_extend_over_the_scrotum_lower_abdomen_or_perineum_answer_more_than_one_if_necessary
    ),
    was_there_a_suspicion_of_urethral_injury_at_presentation = as.factor(was_there_a_suspicion_of_urethral_injury_at_presentation),
    if_there_was_clinical_suspicion_of_a_urethral_injury_why_was_this_select_all_that_apply = as.factor(
      if_there_was_clinical_suspicion_of_a_urethral_injury_why_was_this_select_all_that_apply
    ),
    was_imaging_performed_if_so_what_imaging_select_all_that_apply = as.factor(
      was_imaging_performed_if_so_what_imaging_select_all_that_apply
    ),
    was_the_patient_transferred_between_hospitals_for_specialist_imaging = as.factor(
      was_the_patient_transferred_between_hospitals_for_specialist_imaging
    ),
    did_the_reporting_radiologist_have_a_specialist_interest_in_urology = as.factor(
      did_the_reporting_radiologist_have_a_specialist_interest_in_urology
    ),
    what_time_did_surgical_repair_take_place_knife_to_skin = case_when(
      str_detect(what_time_did_surgical_repair_take_place_knife_to_skin, ":") ~
        dmy_hm(what_time_did_surgical_repair_take_place_knife_to_skin),
      TRUE ~ dmy(what_time_did_surgical_repair_take_place_knife_to_skin) + hours(12)
    ),
    did_the_patient_get_to_theatre_within_24_hours_of_injury = as.factor(did_the_patient_get_to_theatre_within_24_hours_of_injury),
    did_the_patient_get_to_theatre_within_24_hours_of_presentation = as.factor(
      did_the_patient_get_to_theatre_within_24_hours_of_presentation
    ),
    was_the_patient_transferred_from_another_hospital_for_surgical_repair = as.factor(
      was_the_patient_transferred_from_another_hospital_for_surgical_repair
    ),
    were_antibiotics_given_at_induction = as.factor(were_antibiotics_given_at_induction),
    was_a_consultant_present_in_theatre = as.factor(was_a_consultant_present_in_theatre),
    was_a_trainee_present_in_theatre = as.factor(was_a_trainee_present_in_theatre),
    what_was_the_grade_of_the_primary_surgeon = as.factor(what_was_the_grade_of_the_primary_surgeon),
    was_a_specialist_opinion_sought_from_a_urologist_with_a_specialist_interest_in_andrology = as.factor(
      was_a_specialist_opinion_sought_from_a_urologist_with_a_specialist_interest_in_andrology
    ),
    was_a_urologist_with_a_specialist_interest_in_andrology_present_in_theatre = as.factor(
      was_a_urologist_with_a_specialist_interest_in_andrology_present_in_theatre
    ),
    where_was_the_site_of_the_fracture_answer_all_that_apply = as.factor(where_was_the_site_of_the_fracture_answer_all_that_apply),
    what_incision_was_made = as.factor(what_incision_was_made),
    was_a_concurrent_circumcision_performed = as.factor(was_a_concurrent_circumcision_performed),
    what_suture_was_used_for_the_repair_of_the_tunical_fracture = as.factor(
      what_suture_was_used_for_the_repair_of_the_tunical_fracture
    ),
    were_intraoperative_urethral_investigations_performed_select_all_that_apply = as.factor(
      were_intraoperative_urethral_investigations_performed_select_all_that_apply
    ),
    was_a_urethral_injury_noted_at_surgery = as.factor(was_a_urethral_injury_noted_at_surgery),
    if_there_was_a_urethral_injury_what_suture_was_used_to_repair_the_urethra = as.factor(
      if_there_was_a_urethral_injury_what_suture_was_used_to_repair_the_urethra
    ),
    were_there_any_early_complications_clavien_dindo_3 = as.factor(were_there_any_early_complications_clavien_dindo_3),
    if_early_complications_occurred_what_were_they = as.factor(if_early_complications_occurred_what_were_they),
    if_other_complications_occurred_please_specify_enter_na_if_not_relevant = as.factor(
      if_other_complications_occurred_please_specify_enter_na_if_not_relevant
    ),
    was_follow_up_arranged = as.factor(was_follow_up_arranged),
    when_was_follow_up_if_organised_if_no_follow_up_enter_01_01_2000 = dmy(
      when_was_follow_up_if_organised_if_no_follow_up_enter_01_01_2000
    ),
    was_follow_up_arranged_with_a_urologist_with_an_interest_in_andrology = as.factor(
      was_follow_up_arranged_with_a_urologist_with_an_interest_in_andrology
    ),
    presentation_to_knife_to_skin = hms(presentation_to_knife_to_skin),
    
    # Outcomes
    at_follow_up_did_the_patient_describe_ongoing_pain = factor(
      at_follow_up_did_the_patient_describe_ongoing_pain,
      levels = c("Yes", "No", "No Follow up", "Not from the information available")
    ),
    did_the_patient_describe_new_or_worsening_erectile_dysfunction_at_follow_up = factor(
      did_the_patient_describe_new_or_worsening_erectile_dysfunction_at_follow_up,
      levels = c(
        "Yes - Unable to quantify from information available",
        "Yes - Manageable with medication",
        "Yes - Still bothersome despite medication",
        "No",
        "No follow up"
      )
    ),
    at_follow_up_did_the_patient_describe_penile_curvature = factor(
      at_follow_up_did_the_patient_describe_penile_curvature,
      levels = c(
        "Yes - Able to have penetrative intercourse",
        "Yes - Not clear if limiting penetrative intercourse",
        "No",
        "No follow up"
      )
    ),
    at_follow_up_was_there_evidence_of_penile_wasting_or_shortening = factor(
      at_follow_up_was_there_evidence_of_penile_wasting_or_shortening,
      levels = c("Yes", "No", "No follow up")
    ),
    laterality = case_when(
      str_detect(where_was_the_site_of_the_fracture_answer_all_that_apply, regex("left", ignore_case = TRUE)) &
        !str_detect(where_was_the_site_of_the_fracture_answer_all_that_apply, regex("right", ignore_case = TRUE)) ~ "Unilateral",
      str_detect(where_was_the_site_of_the_fracture_answer_all_that_apply, regex("right", ignore_case = TRUE)) &
        !str_detect(where_was_the_site_of_the_fracture_answer_all_that_apply, regex("left", ignore_case = TRUE)) ~ "Unilateral",
      str_detect(where_was_the_site_of_the_fracture_answer_all_that_apply, regex("left", ignore_case = TRUE)) &
        str_detect(where_was_the_site_of_the_fracture_answer_all_that_apply, regex("right", ignore_case = TRUE)) ~ "Bilateral",
      TRUE ~ NA_character_
    ),
    n_locations = sum(c(
      str_detect(where_was_the_site_of_the_fracture_answer_all_that_apply, "Distal"),
      str_detect(where_was_the_site_of_the_fracture_answer_all_that_apply, "Mid"),
      str_detect(where_was_the_site_of_the_fracture_answer_all_that_apply, "Proximal"),
      str_detect(where_was_the_site_of_the_fracture_answer_all_that_apply, "Crura")
    )),
    locality = case_when(
      n_locations >= 2 ~ "Multi-level",
      str_detect(where_was_the_site_of_the_fracture_answer_all_that_apply, "Distal") ~ "Distal",
      str_detect(where_was_the_site_of_the_fracture_answer_all_that_apply, "Mid") ~ "Mid",
      str_detect(where_was_the_site_of_the_fracture_answer_all_that_apply, "Proximal") ~ "Proximal",
      str_detect(where_was_the_site_of_the_fracture_answer_all_that_apply, "Crura") ~ "Crura",
      TRUE ~ NA_character_
    ),
    
    # Other
    are_there_any_other_comments_you_would_like_to_make_regarding_the_case = as.character(
      are_there_any_other_comments_you_would_like_to_make_regarding_the_case
    ),
    .keep = "all"
  ) %>%
  select(
    what_was_the_estimated_date_and_time_of_injury,
    when_was_the_patient_first_seen_in_hospital_by_a_healthcare_professional,
    what_time_did_surgical_repair_take_place_knife_to_skin,
    everything()
  ) %>% 
  select(-n_locations)

## Warning: There were 799 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `what_was_the_estimated_date_and_time_of_injury =
##   case_when(...)`.
## ℹ In row 1.
## Caused by warning:
## ! All formats failed to parse. No formats found.
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 798 remaining warnings.

2.9 Define outcome data

outcomes <- snap %>%
  subset(
    select = c(
      submission_id,
      at_follow_up_did_the_patient_describe_ongoing_pain,
      did_the_patient_describe_new_or_worsening_erectile_dysfunction_at_follow_up,
      at_follow_up_did_the_patient_describe_penile_curvature,
      at_follow_up_was_there_evidence_of_penile_wasting_or_shortening,
      were_there_any_early_complications_clavien_dindo_3
    )
  ) %>% 
  mutate(good_outcome = case_when(
    at_follow_up_did_the_patient_describe_ongoing_pain == "No" & 
      did_the_patient_describe_new_or_worsening_erectile_dysfunction_at_follow_up == "No" & 
      at_follow_up_did_the_patient_describe_penile_curvature == "No" & 
      at_follow_up_was_there_evidence_of_penile_wasting_or_shortening == "No" & 
      were_there_any_early_complications_clavien_dindo_3 == "No" ~ "Yes",
    at_follow_up_did_the_patient_describe_ongoing_pain == "Yes" |
      did_the_patient_describe_new_or_worsening_erectile_dysfunction_at_follow_up == "Yes" | 
      at_follow_up_did_the_patient_describe_penile_curvature == "Yes" |
      at_follow_up_was_there_evidence_of_penile_wasting_or_shortening == "Yes" |
      were_there_any_early_complications_clavien_dindo_3 == "Yes" ~ "No",
    TRUE ~ NA_character_
  ))

3 Explore Predictors - Descriptive statistics only

3.1 Define Predictors

predictors <- snap %>%
               subset(
                 select = -c(
                   at_follow_up_did_the_patient_describe_ongoing_pain,
                   did_the_patient_describe_new_or_worsening_erectile_dysfunction_at_follow_up,
                   at_follow_up_did_the_patient_describe_penile_curvature,
                   at_follow_up_was_there_evidence_of_penile_wasting_or_shortening,
                   if_other_complications_occurred_please_specify_enter_na_if_not_relevant,
                   if_early_complications_occurred_what_were_they,
                   were_there_any_early_complications_clavien_dindo_3,
                   are_there_any_other_comments_you_would_like_to_make_regarding_the_case,
                   was_follow_up_arranged_with_a_urologist_with_an_interest_in_andrology,
                   when_was_follow_up_if_organised_if_no_follow_up_enter_01_01_2000
                 )
               )

3.2 Plot missing

plot_missing(predictors)

3.3 Age

3.3.1 Histogram

ggplot(predictors, aes(x = what_age_was_the_patient_at_the_time_of_the_injury)) +
  geom_histogram(bins = 50, fill = "skyblue", color = "black") +
  labs(
    title = "Histogram of Age at Time of Injury",
    x = "Age (years)",
    y = "Count"
  ) +
  theme_minimal()

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

3.3.2 Density plot

ggplot(predictors, aes(x = what_age_was_the_patient_at_the_time_of_the_injury)) +
  geom_density(fill = "skyblue", color = "black") +
  labs(
    title = "Density Plot of Age at Time of Injury",
    x = "Age (years)",
    y = "Count"
  ) +
  theme_minimal()

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_density()`).

3.4 Type of Suture

3.4.1 Overall

predictors %>% 
  drop_na(what_suture_was_used_for_the_repair_of_the_tunical_fracture) %>%
  tabyl(what_suture_was_used_for_the_repair_of_the_tunical_fracture) %>%
  mutate(
    "Type of Suture" = what_suture_was_used_for_the_repair_of_the_tunical_fracture,
    Percent = round(percent * 100),
         .keep = "unused") %>%
  select(c("Type of Suture",
             n,
             Percent)) %>%
  gt() %>%
  gt_theme_espn()

Type of Suture	n	Percent
PDS 3-0 Suture	46	19
PDS 2-0 Suture	93	38
Vicryl 3-0 Suture	5	2
Other	30	12
Nylon 4-0 Suture	1	0
PDS 0-Suture	37	15
Vicryl 4-0 Suture	3	1
PDS 1-0 Suture	4	2
PDS 4-0 Suture	10	4
Vicryl 2-0 Suture	10	4
Vicryl 0 -Suture	2	1
PDS 2-0 Suture - Vicryl 3-0 Suture	1	0
PDS 0-Suture - PDS 2-0 Suture	2	1
Nylon 3-0 Suture	1	0
PDS 2-0 Suture - Vicryl 2-0 Suture	1	0

predictors <- predictors %>%
  mutate(
    type_of_suture = case_when(
      str_detect(
        what_suture_was_used_for_the_repair_of_the_tunical_fracture,
        "PDS"
      ) ~ "Non-absorbable",
      str_detect(
        what_suture_was_used_for_the_repair_of_the_tunical_fracture,
        "Nylon"
      ) ~ "Non-absorbable",
      str_detect(
        what_suture_was_used_for_the_repair_of_the_tunical_fracture,
        "Vicryl"
      ) ~ "Absorbable",
      what_suture_was_used_for_the_repair_of_the_tunical_fracture == "Other" |
        what_suture_was_used_for_the_repair_of_the_tunical_fracture == "PDS 2-0 Suture - Vicryl 2-0 Suture" ~ "Unclear",
      TRUE ~ NA_character_
    )
  )

3.4.2 Simplified into absorbable vs non-absorbable

predictors %>% 
  select(type_of_suture) %>%
  group_by(type_of_suture) %>%
  summarise(
    n = n(),
    "Percentage (%)" = round(n / 258 * 100)
  ) %>%
  gt() %>%
  gt_theme_espn()

type_of_suture	n	Percentage (%)
Absorbable	20	8
Non-absorbable	196	76
Unclear	30	12
NA	12	5

3.5 Urethral Injury

predictors %>% 
  drop_na(was_a_urethral_injury_noted_at_surgery) %>%
  tabyl(was_a_urethral_injury_noted_at_surgery) %>%
  mutate(
    "Urethral Injury?" = was_a_urethral_injury_noted_at_surgery,
    Percent = round(percent * 100),
         .keep = "unused") %>%
  select(c("Urethral Injury?",
             n,
             Percent)) %>%
  gt() %>%
  gt_theme_espn()

Urethral Injury?	n	Percent
No	200	80
Yes	51	20

3.6 History of Prior Penile Fracture

predictors %>% 
  drop_na(did_the_patient_have_a_history_of_a_prior_penile_fracture) %>%
  tabyl(did_the_patient_have_a_history_of_a_prior_penile_fracture) %>%
  mutate(
    "Prior Penile Fracture?" = did_the_patient_have_a_history_of_a_prior_penile_fracture,
    Percent = round(percent * 100),
         .keep = "unused") %>%
  select(c("Prior Penile Fracture?",
             n,
             Percent)) %>%
  gt() %>%
  gt_theme_espn()

Prior Penile Fracture?	n	Percent
No	233	93
Yes	18	7

3.7 Mechanism of Injury

predictors %>% 
  drop_na(mechanism_of_injury) %>%
  tabyl(mechanism_of_injury) %>%
  mutate(
    "Mechanism of Injury" = mechanism_of_injury,
    Percent = round(percent * 100),
         .keep = "unused") %>%
  select(c("Mechanism of Injury",
             n,
             Percent)) %>%
  gt() %>%
  gt_theme_espn()

Mechanism of Injury	n	Percent
Sexual intercourse trauma	214	86
Other	6	2
Non-sexual blunt force trauma	19	8
Masturbation injury	10	4

3.8 If Sexual Intercourse - What position?

This could be amalgamated with Mechanism of injury

predictors %>% 
  drop_na(if_the_injury_was_sustained_during_sexual_intercourse_what_sexual_position_was_involved) %>%
  filter(if_the_injury_was_sustained_during_sexual_intercourse_what_sexual_position_was_involved != "NA - Injury not sustained during sexual intercourse") %>%
  mutate(if_the_injury_was_sustained_during_sexual_intercourse_what_sexual_position_was_involved = factor(if_the_injury_was_sustained_during_sexual_intercourse_what_sexual_position_was_involved,
                                                                                                          levels = c(
                                                                                                            "Patient on top (Missionary/face to face)",
                                                                                                            "Partner on top of patient",
                                                                                                            "Partner on top of patient and facing away (Reverse cowgirl/cowboy)",
                                                                                                            "Patient behind partner (Doggy style)",
                                                                                                            "Side to side (spooning)",
                                                                                                            "Other"
                                                                                                          ))) %>%
  tabyl(if_the_injury_was_sustained_during_sexual_intercourse_what_sexual_position_was_involved) %>%
  mutate(
    "Position" = if_the_injury_was_sustained_during_sexual_intercourse_what_sexual_position_was_involved,
    Percent = round(percent * 100),
         .keep = "unused") %>%
  select(c("Position",
             n,
             Percent)) %>%
  gt() %>%
  gt_theme_espn()

Position	n	Percent
Patient on top (Missionary/face to face)	17	30
Partner on top of patient	13	23
Partner on top of patient and facing away (Reverse cowgirl/cowboy)	3	5
Patient behind partner (Doggy style)	21	38
Side to side (spooning)	1	2
Other	1	2

3.9 Pain at Presentation

predictors %>% 
  drop_na(was_pain_present_at_presentation) %>%
  tabyl(was_pain_present_at_presentation) %>%
  mutate(
    "Pain at Presentation" = was_pain_present_at_presentation,
    Percent = round(percent * 100),
         .keep = "unused") %>%
  select(c("Pain at Presentation",
             n,
             Percent)) %>%
  gt() %>%
  gt_theme_espn()

Pain at Presentation	n	Percent
No	35	15
Yes	191	85

3.10 Detumescence

predictors %>% 
  drop_na(was_there_immediate_detumescence_loss_of_erection) %>%
  tabyl(was_there_immediate_detumescence_loss_of_erection) %>%
  mutate(
    "Detumescence" = was_there_immediate_detumescence_loss_of_erection,
    Percent = round(percent * 100),
         .keep = "unused") %>%
  select(c("Detumescence",
             n,
             Percent)) %>%
  gt() %>%
  gt_theme_espn()

Detumescence	n	Percent
Yes	179	88
No	25	12

3.11 Audible Pop

predictors %>% 
  drop_na(was_an_audible_pop_reported_by_the_patient) %>%
  tabyl(was_an_audible_pop_reported_by_the_patient) %>%
  mutate(
    "Audible Pop" = was_an_audible_pop_reported_by_the_patient,
    Percent = round(percent * 100),
         .keep = "unused") %>%
  select(c("Audible Pop",
             n,
             Percent)) %>%
  gt() %>%
  gt_theme_espn()

Audible Pop	n	Percent
No	61	32
Yes	130	68

3.12 Bruising

Composite variable based on ‘was_penile_bruising_present’ and  ‘did_bruising_extend_over_the_scrotum_lower_abdomen_or_perineum_answer_more_than_one_if_necessary’

bruise_summary <- predictors %>%
  drop_na(was_penile_bruising_present) %>%
  mutate(
    scrotum = case_when(
      did_bruising_extend_over_the_scrotum_lower_abdomen_or_perineum_answer_more_than_one_if_necessary %in% c(
        "Scrotum",
        "Scrotum  -  Lower abdomen",
        "Scrotum  -  Lower abdomen  -  Perineum",
        "Scrotum  -  Not clear from information available",
        "Scrotum  -  Perineum",
        "Scrotum  -  Lower abdomen  -  Perineum  -  No bruising present  -  Not clear from information available"
      ) ~ "Yes",
      did_bruising_extend_over_the_scrotum_lower_abdomen_or_perineum_answer_more_than_one_if_necessary == "No bruising present" |
        was_penile_bruising_present == "No" ~ "No",
      TRUE ~ NA_character_
    ),
    lower_abdomen = case_when(
      did_bruising_extend_over_the_scrotum_lower_abdomen_or_perineum_answer_more_than_one_if_necessary %in% c(
        "Lower abdomen",
        "Scrotum  -  Lower abdomen",
        "Scrotum  -  Lower abdomen  -  Perineum",
        "Scrotum  -  Lower abdomen  -  Perineum  -  No bruising present  -  Not clear from information available"
      ) ~ "Yes",
      did_bruising_extend_over_the_scrotum_lower_abdomen_or_perineum_answer_more_than_one_if_necessary == "No bruising present" |
        was_penile_bruising_present == "No" ~ "No",
      TRUE ~ NA_character_
    ),
    perineum = case_when(
      did_bruising_extend_over_the_scrotum_lower_abdomen_or_perineum_answer_more_than_one_if_necessary %in% c(
        "Perineum",
        "Scrotum  -  Lower abdomen  -  Perineum",
        "Scrotum  -  Perineum",
        "Scrotum  -  Lower abdomen  -  Perineum  -  No bruising present  -  Not clear from information available"
      ) ~ "Yes",
      did_bruising_extend_over_the_scrotum_lower_abdomen_or_perineum_answer_more_than_one_if_necessary == "No bruising present" |
        was_penile_bruising_present == "No" ~ "No",
      TRUE ~ NA_character_
    ),
    no_bruising = case_when(
      scrotum == "No" & lower_abdomen == "No" & perineum == "No" ~ "Yes",
      scrotum == "Yes" | lower_abdomen == "Yes" | perineum == "Yes" ~ "No",
      TRUE ~ NA_character_
    )
  ) %>%
  select(scrotum, lower_abdomen, perineum, no_bruising) %>%
  pivot_longer(cols = everything(), names_to = "Location", values_to = "Value") %>%
  group_by(Location, Value) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(Location) %>%
  mutate(pct = 100 * n / sum(n)) %>%
  ungroup() %>%
  mutate(display = paste0(n, " (", sprintf("%.1f", pct), "%)")) %>%
  select(Location, Value, display) %>%
  pivot_wider(names_from = Value, values_from = display)


bruise_summary %>%
  kable(
    caption = "Counts and Percentages of Bruising by Location",
    align = c("l", "c", "c", "c"),  
    col.names = c("Location", "Yes", "No", "Missing")
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = FALSE,
    position = "center"
  ) %>%
  row_spec(0, bold = TRUE, background = "#f2f2f2") %>%  
  column_spec(1, bold = TRUE)

Counts and Percentages of Bruising by Location

Location	Yes	No	Missing
lower_abdomen	106 (42.7%)	21 (8.5%)	121 (48.8%)
no_bruising	89 (35.9%)	106 (42.7%)	53 (21.4%)
perineum	106 (42.7%)	8 (3.2%)	134 (54.0%)
scrotum	106 (42.7%)	79 (31.9%)	63 (25.4%)

3.13 Incision Type

predictors %>% 
  drop_na(what_incision_was_made) %>%
  tabyl(what_incision_was_made) %>%
  mutate(
    "Incision" = what_incision_was_made,
    Percent = round(percent * 100),
         .keep = "unused") %>%
  select(c("Incision",
             n,
             Percent)) %>%
  gt() %>%
  gt_theme_espn()

Incision	n	Percent
longitudinal	59	24
Penoscrotal	96	38
Subcoronal with degloving of penis	81	32
Subcoronal without degloving of penis	7	3
Other	8	3

3.14 Transfer from Another Hospital

predictors %>% 
  drop_na(was_the_patient_transferred_from_another_hospital_for_surgical_repair) %>%
  tabyl(was_the_patient_transferred_from_another_hospital_for_surgical_repair) %>%
  mutate(
    "Transfer for Repair" = was_the_patient_transferred_from_another_hospital_for_surgical_repair,
    Percent = round(percent * 100),
         .keep = "unused") %>%
  select(c("Transfer for Repair",
             n,
             Percent)) %>%
  gt() %>%
  gt_theme_espn()

Transfer for Repair	n	Percent
No	198	77
Yes	58	23

3.15 Time Variables

3.15.1 Time from Presentation to Theatre

predictors <- predictors %>% mutate(
  presentation_to_knife_to_skin = difftime(
    what_time_did_surgical_repair_take_place_knife_to_skin,
    when_was_the_patient_first_seen_in_hospital_by_a_healthcare_professional,
    units = "mins"
  ) %>% as.numeric()
)


presentation_to_knife_to_skin_plot <- ggplot(predictors, aes(x = presentation_to_knife_to_skin / 60)) +
  geom_histogram(bins = 50, fill = "skyblue", color = "black") +
  labs(
    title = "Histogram of Time from Presentation to Theatre",
    x = "Hours",
    y = "Count"
  ) + xlim(0,50) +
  theme_minimal() 
  
presentation_to_knife_to_skin_plot  

## Warning: Removed 56 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_bar()`).

3.15.2 Time from Time of Injury to Theatre

time_from_injury_to_th <- predictors %>% 
  mutate(
  time_from_injury_to_th = difftime(what_time_did_surgical_repair_take_place_knife_to_skin,
                                    what_was_the_estimated_date_and_time_of_injury,
                                  units = "mins")
) %>%
  mutate(
  time_from_injury_to_th = ifelse(time_from_injury_to_th < 0 | time_from_injury_to_th > (140*60),
                                  NA,
                                  time_from_injury_to_th)
)

time_from_injury_to_th$time_from_injury_to_th <- as.numeric(time_from_injury_to_th$time_from_injury_to_th)

ggplot(time_from_injury_to_th, aes(x = time_from_injury_to_th / 60)) +
  geom_histogram(bins = 50, fill = "skyblue", color = "black") +
  labs(
    title = "Histogram of Time from Injury to Theatre",
    x = "Hours",
    y = "Count"
  ) +
  theme_minimal()

## Warning: Removed 41 rows containing non-finite outside the scale range
## (`stat_bin()`).

3.15.3 Transform Time variables into Bins

time_from_injury_to_th <- time_from_injury_to_th %>%
  mutate(
    time_from_injury_to_th1 = time_from_injury_to_th / 60,
    presentation_to_knife_to_skin1 = presentation_to_knife_to_skin / 60,
    time_from_injury_to_th_bins = cut(
      time_from_injury_to_th1,
      breaks = c(-Inf, 6, 12, 24, 36, Inf),
      labels = c("<6h", "6–12h", "12–24h", "24-36h", ">36h")
    ),
    presentation_to_knife_to_skin_bins = cut(
      presentation_to_knife_to_skin1,
      breaks = c(-Inf, 6, 12, 24, 36, Inf),
      labels = c("<6h", "6–12h", "12–24h", "24-36h", ">36h")
    )
  )

p1 <- time_from_injury_to_th %>% 
  drop_na(time_from_injury_to_th_bins) %>%
  ggplot(aes(x = time_from_injury_to_th_bins)) +
  geom_bar(fill = "#2C7BB6", colour = "black") +
  labs(
    x = "Time from Injury to Theatre (hours, binned)",
    y = "Count"
  ) + ylim(0,100) +
  theme_minimal()

p2 <- time_from_injury_to_th %>% 
  drop_na(presentation_to_knife_to_skin_bins) %>%
  ggplot(aes(x = presentation_to_knife_to_skin_bins)) +
  geom_bar(fill = "#D7191C", colour = "black") +
  labs(
    x = "Presentation to Knife-to-Skin (hours, binned)",
    y = "Count"
  ) + ylim(0,100) +
  theme_minimal()

time_binned_plot <- plot_grid(p1, p2, labels = c("A", "B"))
time_binned_plot

3.16 Imaging

3.16.1 Correlation between Imaging Performed and Reporting Radiologist Subspecialism

3.16.1.1 Any imaging

predictors <- predictors %>%
  mutate(
    was_imaging_performed_if_so_what_imaging_select_all_that_apply = recode_factor(was_imaging_performed_if_so_what_imaging_select_all_that_apply,
      "No" = "No",
      "Yes - Both Ultrasound and MRI  -  Not clear from information available" = NA_character_,
      "Yes - Both Ultrasound and MRI" = "MRI & US",
      "Yes - MRI penis" = "MRI",
      "Yes - MRI penis  -  Yes - Ultrasound" = "MRI & US",
      "Yes - Ultrasound" = "US"
    ),
    did_the_reporting_radiologist_have_a_specialist_interest_in_urology = recode_factor(did_the_reporting_radiologist_have_a_specialist_interest_in_urology,
      "No" = "No",
      "Yes" = "Yes",
      "Not clear from information available" = NA_character_
    ),
    was_the_patient_transferred_between_hospitals_for_specialist_imaging = recode_factor(was_the_patient_transferred_between_hospitals_for_specialist_imaging,
      "No" = "No",
      "Yes" = "Yes",
      "Not clear from information available" = NA_character_
    )
  )

tab <- predictors %>%
  drop_na(was_imaging_performed_if_so_what_imaging_select_all_that_apply,
          did_the_reporting_radiologist_have_a_specialist_interest_in_urology) %>%
  mutate(any_imaging = ifelse(was_imaging_performed_if_so_what_imaging_select_all_that_apply == "No",
                              "No",
                              "Yes") %>% as.factor()) %>%
  tabyl(
    any_imaging,
    did_the_reporting_radiologist_have_a_specialist_interest_in_urology
  ) %>%
  adorn_totals("row") %>%            
  adorn_totals("col") %>%            
  adorn_percentages("row") 

# Get counts separately
tab_counts <- predictors %>%
  drop_na(was_imaging_performed_if_so_what_imaging_select_all_that_apply,
          did_the_reporting_radiologist_have_a_specialist_interest_in_urology) %>%
  tabyl(was_imaging_performed_if_so_what_imaging_select_all_that_apply,
          did_the_reporting_radiologist_have_a_specialist_interest_in_urology) %>%
  adorn_totals("row") %>%
  adorn_totals("col")

# Combine counts and percentages into n (%)
tab_n_pct <- tab_counts
for(col in 2:ncol(tab_counts)) {
  tab_n_pct[[col]] <- paste0(tab_counts[[col]], " (", 
                             round(tab[[col]]*100, 1), "%)")
}

tab_n_pct %>%
  kable("html", escape = FALSE, 
        col.names = c("Imaging Type", "Not Reported by Subspecialist Radiologist", "Reported by Subspecialist Radiologist", "Total")) %>%
  kable_styling(full_width = FALSE, 
                bootstrap_options = c("striped", "hover", "condensed")) %>%
  row_spec(0, bold = TRUE, background = "#D3D3D3") 

Imaging Type	Not Reported by Subspecialist Radiologist	Reported by Subspecialist Radiologist	Total
No	52 (98.1%)	1 (1.9%)	53 (100%)
MRI & US	0 (14%)	10 (86%)	10 (100%)
MRI	11 (39.7%)	36 (60.3%)	47 (100%)
US	6 (98.1%)	58 (1.9%)	64 (100%)
Total	69 (14%)	105 (86%)	174 (100%)

3.16.1.2 Specific Imaging

tab <- predictors %>%
  drop_na(was_imaging_performed_if_so_what_imaging_select_all_that_apply,
          did_the_reporting_radiologist_have_a_specialist_interest_in_urology) %>%
  tabyl(
    was_imaging_performed_if_so_what_imaging_select_all_that_apply,
    did_the_reporting_radiologist_have_a_specialist_interest_in_urology
  ) %>%
  adorn_totals("row") %>%            
  adorn_totals("col") %>%            
  adorn_percentages("row") 

# Get counts separately
tab_counts <- predictors %>%
  drop_na(was_imaging_performed_if_so_what_imaging_select_all_that_apply,
          did_the_reporting_radiologist_have_a_specialist_interest_in_urology) %>%
  tabyl(was_imaging_performed_if_so_what_imaging_select_all_that_apply,
          did_the_reporting_radiologist_have_a_specialist_interest_in_urology) %>%
  adorn_totals("row") %>%
  adorn_totals("col")

# Combine counts and percentages into n (%)
tab_n_pct <- tab_counts
for(col in 2:ncol(tab_counts)) {
  tab_n_pct[[col]] <- paste0(tab_counts[[col]], " (", 
                             round(tab[[col]]*100, 1), "%)")
}

tab_n_pct %>%
  kable("html", escape = FALSE, 
        col.names = c("Imaging Type", "Not Reported by Subspecialist Radiologist", "Reported by Subspecialist Radiologist", "Total")) %>%
  kable_styling(full_width = FALSE, 
                bootstrap_options = c("striped", "hover", "condensed")) %>%
  row_spec(0, bold = TRUE, background = "#D3D3D3") 

Imaging Type	Not Reported by Subspecialist Radiologist	Reported by Subspecialist Radiologist	Total
No	52 (98.1%)	1 (1.9%)	53 (100%)
MRI & US	0 (0%)	10 (100%)	10 (100%)
MRI	11 (23.4%)	36 (76.6%)	47 (100%)
US	6 (9.4%)	58 (90.6%)	64 (100%)
Total	69 (39.7%)	105 (60.3%)	174 (100%)

3.16.2 Correlation between Imaging Performed and Transfer for Imaging

3.16.2.1 Any Imaging

predictors <- predictors %>% 
  mutate(
    any_imaging = ifelse(was_imaging_performed_if_so_what_imaging_select_all_that_apply == "No",
                              "No",
                              "Yes") %>% as.factor(),
    .keep = "all"
  )


tab <- predictors %>%
  drop_na(was_imaging_performed_if_so_what_imaging_select_all_that_apply,
          was_the_patient_transferred_between_hospitals_for_specialist_imaging) %>%
  mutate(any_imaging = ifelse(was_imaging_performed_if_so_what_imaging_select_all_that_apply == "No",
                              "No",
                              "Yes") %>% as.factor()) %>%
  tabyl(
    any_imaging,
    was_the_patient_transferred_between_hospitals_for_specialist_imaging
  ) %>%
  adorn_totals("row") %>%            
  adorn_totals("col") %>%            
  adorn_percentages("row")         

# Get counts separately
tab_counts <- predictors %>%
  drop_na(was_imaging_performed_if_so_what_imaging_select_all_that_apply,
          was_the_patient_transferred_between_hospitals_for_specialist_imaging) %>%
  tabyl(was_imaging_performed_if_so_what_imaging_select_all_that_apply,
          was_the_patient_transferred_between_hospitals_for_specialist_imaging) %>%
  adorn_totals("row") %>%
  adorn_totals("col")

# Combine counts and percentages into n (%)
tab_n_pct <- tab_counts
for(col in 2:ncol(tab_counts)) {
  tab_n_pct[[col]] <- paste0(tab_counts[[col]], " (", 
                             round(tab[[col]]*100, 1), "%)")
}


tab_n_pct %>%
  kable("html", escape = FALSE, 
        col.names = c("Imaging Type", "No Transfer", "Transfer", "Total")) %>%
  kable_styling(full_width = FALSE, 
                bootstrap_options = c("striped", "hover", "condensed")) %>%
  row_spec(0, bold = TRUE, background = "#D3D3D3") 

Imaging Type	No Transfer	Transfer	Total
No	95 (100%)	0 (0%)	95 (100%)
MRI & US	13 (95.4%)	0 (4.6%)	13 (100%)
MRI	55 (97.2%)	5 (2.8%)	60 (100%)
US	77 (100%)	2 (0%)	79 (100%)
Total	240 (95.4%)	7 (4.6%)	247 (100%)

3.16.2.2 Specific Imaging

tab <- predictors %>%
  drop_na(was_imaging_performed_if_so_what_imaging_select_all_that_apply,
          was_the_patient_transferred_between_hospitals_for_specialist_imaging) %>%
  tabyl(
    was_imaging_performed_if_so_what_imaging_select_all_that_apply,
    was_the_patient_transferred_between_hospitals_for_specialist_imaging
  ) %>%
  adorn_totals("row") %>%            
  adorn_totals("col") %>%            
  adorn_percentages("row") 

# Get counts separately
tab_counts <- predictors %>%
  drop_na(was_imaging_performed_if_so_what_imaging_select_all_that_apply,
          was_the_patient_transferred_between_hospitals_for_specialist_imaging) %>%
  tabyl(was_imaging_performed_if_so_what_imaging_select_all_that_apply,
          was_the_patient_transferred_between_hospitals_for_specialist_imaging) %>%
  adorn_totals("row") %>%
  adorn_totals("col")

# Combine counts and percentages into n (%)
tab_n_pct <- tab_counts
for(col in 2:ncol(tab_counts)) {
  tab_n_pct[[col]] <- paste0(tab_counts[[col]], " (", 
                             round(tab[[col]]*100, 1), "%)")
}

tab_n_pct %>%
  kable("html", escape = FALSE, 
        col.names = c("Imaging Type", "No Transfer", "Transfer", "Total")) %>%
  kable_styling(full_width = FALSE, 
                bootstrap_options = c("striped", "hover", "condensed")) %>%
  row_spec(0, bold = TRUE, background = "#D3D3D3") 

Imaging Type	No Transfer	Transfer	Total
No	95 (100%)	0 (0%)	95 (100%)
MRI & US	13 (100%)	0 (0%)	13 (100%)
MRI	55 (91.7%)	5 (8.3%)	60 (100%)
US	77 (97.5%)	2 (2.5%)	79 (100%)
Total	240 (97.2%)	7 (2.8%)	247 (100%)

3.16.3 Does Getting Imaging correlate with increased time to theatre?

time_from_injury_to_th %>% 
  mutate(
    imaging = case_when(
      was_imaging_performed_if_so_what_imaging_select_all_that_apply == "No" ~ "No",
      str_detect(was_imaging_performed_if_so_what_imaging_select_all_that_apply, "Yes") ~ "Yes",
      TRUE ~ NA_character_
    )
  ) %>%
  select(imaging,
         presentation_to_knife_to_skin,
         was_the_patient_transferred_from_another_hospital_for_surgical_repair) %>%
  drop_na() %>%
  group_by(imaging,
           was_the_patient_transferred_from_another_hospital_for_surgical_repair) %>%
  summarise(
    mean_p2k = mean(presentation_to_knife_to_skin, na.rm = TRUE),
    sd_p2k   = sd(presentation_to_knife_to_skin, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  ggplot(aes(
    x = imaging, 
    y = mean_p2k,
    fill = was_the_patient_transferred_from_another_hospital_for_surgical_repair
  )) +
  geom_col(position = position_dodge()) 

#  geom_errorbar(aes(ymin = mean_p2k - sd_p2k, ymax = mean_p2k + sd_p2k),
#                position = position_dodge(0.9),
#                width = 0.2)

3.16.4 Imaging and Rate of Transfer for Surgical Repair

I wonder if this is type II statistical error

plot_data <- time_from_injury_to_th %>%
  mutate(
    imaging = case_when(
      was_imaging_performed_if_so_what_imaging_select_all_that_apply == "No" ~ "No",
      str_detect(was_imaging_performed_if_so_what_imaging_select_all_that_apply, "Yes") ~ "Yes",
      TRUE ~ NA_character_
    )
  ) %>%
  select(imaging, was_the_patient_transferred_from_another_hospital_for_surgical_repair) %>%
  drop_na() %>%
  count(imaging, was_the_patient_transferred_from_another_hospital_for_surgical_repair) %>%
  group_by(imaging) %>%
  mutate(pct = n / sum(n) * 100)

chi_test <- chisq.test(
  xtabs(n ~ imaging + was_the_patient_transferred_from_another_hospital_for_surgical_repair, data = plot_data)
)

pval_label <- paste0("Chi-square p = ", signif(chi_test$p.value, 3))


# Plot raw counts
p_counts <- ggplot(plot_data, aes(x = imaging, y = n,
                                  fill = was_the_patient_transferred_from_another_hospital_for_surgical_repair)) +
  geom_col(position = position_dodge(width = 0.8), width = 0.7) +
  geom_text(aes(label = n),
            position = position_dodge(width = 0.8), vjust = -0.3, size = 4) +
  annotate("text", x = 1.5, y = max(plot_data$n) * 1.1, label = pval_label,
           size = 5, fontface = "italic") +
  scale_fill_brewer(palette = "Set2", name = "Transferred from another hospital") +
  labs(title = "Raw Counts", x = "Imaging performed", y = "Number of patients") +
  theme_minimal(base_size = 14) +
  theme(plot.title = element_text(face = "bold", hjust = 0.5),
        legend.position = "top")

# Plot percentages
p_pct <- ggplot(plot_data, aes(x = imaging, y = pct,
                               fill = was_the_patient_transferred_from_another_hospital_for_surgical_repair)) +
  geom_col(position = position_dodge(width = 0.8), width = 0.7) +
  geom_text(aes(label = paste0(round(pct, 1), "%")),
            position = position_dodge(width = 0.8), vjust = -0.3, size = 4) +
  scale_fill_brewer(palette = "Set2", name = "Transferred from another hospital") +
  labs(title = "Percentages", x = "Imaging performed", y = "Percentage of patients") +
  theme_minimal(base_size = 14) +
  theme(plot.title = element_text(face = "bold", hjust = 0.5),
        legend.position = "top")

plot_grid(p_counts, p_pct)

3.16.5 Imaging and Circumcision Rate

plot_data <- time_from_injury_to_th %>%
  mutate(
    imaging = case_when(
      was_imaging_performed_if_so_what_imaging_select_all_that_apply == "No" ~ "No",
      str_detect(was_imaging_performed_if_so_what_imaging_select_all_that_apply, "Yes") ~ "Yes",
      TRUE ~ NA_character_
    )
  ) %>%
  select(imaging, was_a_concurrent_circumcision_performed) %>%
  filter(was_a_concurrent_circumcision_performed != "Patient previously circumcised") %>%
  drop_na() %>%
  count(imaging, was_a_concurrent_circumcision_performed) %>%
  group_by(imaging) %>%
  mutate(pct = n / sum(n) * 100)

plot_data$was_a_concurrent_circumcision_performed <- factor(plot_data$was_a_concurrent_circumcision_performed,
                                                            levels = c("Yes",
                                                                       "No"))

chi_test <- chisq.test(
  xtabs(n ~ imaging + was_a_concurrent_circumcision_performed, data = plot_data)
)

pval_label <- paste0("Chi-square p = ", signif(chi_test$p.value, 3))


# Plot raw counts
p_counts <- ggplot(plot_data, aes(x = imaging, y = n,
                                  fill = was_a_concurrent_circumcision_performed)) +
  geom_col(position = position_dodge(width = 0.8), width = 0.7) +
  geom_text(aes(label = n),
            position = position_dodge(width = 0.8), vjust = -0.3, size = 4) +
  annotate("text", x = 1.5, y = max(plot_data$n) * 1.1, label = pval_label,
           size = 5, fontface = "italic") +
  scale_fill_brewer(palette = "Set2", name = "Circumcision Performed") +
  labs(title = "Raw Counts", x = "Imaging performed", y = "Number of patients") +
  theme_minimal(base_size = 14) +
  theme(plot.title = element_text(face = "bold", hjust = 0.5),
        legend.position = "top")

# Plot percentages
p_pct <- ggplot(plot_data, aes(x = imaging, y = pct,
                               fill = was_a_concurrent_circumcision_performed)) +
  geom_col(position = position_dodge(width = 0.8), width = 0.7) +
  geom_text(aes(label = paste0(round(pct, 1), "%")),
            position = position_dodge(width = 0.8), vjust = -0.3, size = 4) +
  scale_fill_brewer(palette = "Set2", name = "Circumcision Performed") +
  labs(title = "Percentages", x = "Imaging performed", y = "Percentage of patients") +
  theme_minimal(base_size = 14) +
  theme(plot.title = element_text(face = "bold", hjust = 0.5),
        legend.position = "top")

plot_grid(p_counts, p_pct)

3.17 Volume

3.17.1 How many repairs have individual units undertaken?

volume <- snap %>%
  select(submission_id,
         please_select_the_deanery_hospital_you_are_completing_the_form_for) %>%
  group_by(please_select_the_deanery_hospital_you_are_completing_the_form_for) %>%
  summarise(n = n())

volume %>% 
  arrange(desc(n)) %>%
  gt() %>% gt_theme_espn()

please_select_the_deanery_hospital_you_are_completing_the_form_for	n
London, North University College Hospital, London	39
London, South King's College Hospital, London	10
Yorkshire & Humber St James University Hospital, Leeds	10
East Midlands Nottingham City Hospital	8
South West Southmead Hospital, Bristol	8
North East Freeman Hospital, Newcastle	7
London, South St George's Hospital, London	6
Scotland Western General Hospital, Edinburgh	6
South Central Southampton General Hospital	6
East Midlands Leicester General Hospital	5
East of England Norfolk and Norwich University Hospital	5
London, South Guy's and Thomas's Hospital	5
South Central Queen Alexandra Hospital, Portsmouth	5
South West Cheltenham General Hospital	5
South West Royal Devon & Exeter Hospital	5
East Midlands Royal Derby Hospital	4
Kent, Surrey, Sussex Princess Royal Hospital, Hayward's Heath	4
London, North Northwick Park Hospital	4
North East Sunderland Royal Hospital	4
North West Cumberland Infirmary, Carlisle	4
Yorkshire & Humber Castle Hill Hospital, Hull	4
Kent, Surrey, Sussex Frimley Park Hospital	3
London, North Barking, Havering & Redbridge University Hospital	3
London, North West Middlesex University Hospital	3
London, North Whipps Cross Hospital, London	3
North West Arrowe Parke Hospital, Wirral	3
North West Blackpool Victoria Hospital	3
North West Royal Liverpool & Broadgreen University Hospitals	3
North West Stepping Hill Hospital, Stockport	3
North West Whiston Hospital	3
South Central Royal Berkshire Hospital	3
South West Derriford Hospital, Plymouth	3
South West Royal Cornwall Hospital, Truro	3
South West The Great Western Hospital, Swindon	3
Wales University Hospital of Wales, Cardiff	3
West Midlands Royal Stoke University Hospitals	3
West Midlands Sandwell District General Hospital, West Bromwich	3
Yorkshire & Humber Bradford Royal Infirmary	3
Yorkshire & Humber Pinderfields General Hospital, Wakefield	3
East of England Ipswich Hospital, Suffolk	2
East of England Southend Hospital	2
Kent, Surrey, Sussex Kent & Canterbury Hospital	2
Kent, Surrey, Sussex Medway Maritime Hospital, Gillingham	2
North East The James Cook University Hospital, Middlesborough	2
North West Manchester Royal Infirmary	2
North West North Manchester General Hospital	2
North West Royal Blackburn Hospital	2
North West Royal Bolton Hospital	2
North West Warrington District General Hospital	2
Scotland Glasgow Royal Infirmary	2
South West Musgrove Park Hospital, Taunton	2
Wales Royal Gwent Hospital	2
West Midlands Alexandra Hospital, Redditch	2
West Midlands University Hospital, Coventry	2
East Midlands Chesterfield Royal Hospital	1
East Midlands Kettering General Hospital	1
East Midlands Lincoln County Hospital	1
East Midlands Northampton General Hospital	1
East of England Addenbrooke's Hospital, Cambridge	1
East of England Peterborough City Hospital	1
East of England Queen Elizabeth Hospital, King's Lynn	1
Kent, Surrey, Sussex Eastbourne District General Hospital	1
London, North Chelsea & Westminster Hospital, London	1
London, North Hillingdon Hospital, Uxbridge	1
North East Bishop Auckland Hospital	1
Scotland Aberdeen Royal Infirmary	1
Scotland Queen Elizabeth University Hospital, Glasgow	1
South West Royal United Hospital, Bath	1
Wales Morriston Hospital, Swansea	1
West Midlands Hereford County Hospital	1
West Midlands Queen's Hospital, Burton-on-Trent	1
Yorkshire & Humber Airedale Hospital	1
Yorkshire & Humber Royal Hallamshire Hospital, Sheffield	1
Yorkshire & Humber Scunthorpe General Hospital	1

3.17.2 Subdivide into Low/Medium/High Volume centres

UCLH seems to be an outlier in terms of volume - this should be “High”
Then arbitrarily subdivide into “Medium” (>5) and “Low” (</=5)

volume <- volume %>% 
  ungroup() %>%
  mutate(
    volume = case_when(
      n > 30 ~ "High",
      n > 5 & n <= 30 ~ "Medium",
      n < 6 ~ "Low" 
    ),
    volume_n = n
  )

volume$volume <- factor(volume$volume,
                        levels = c("High",
                                   "Medium",
                                   "Low"))

volume_clean <- volume %>%
  mutate(hospital_clean = str_squish(str_replace_all(
    please_select_the_deanery_hospital_you_are_completing_the_form_for, "[\r\n]", " "
  )))

predictors <- predictors %>% 
  mutate(hospital_clean = str_squish(str_replace_all(
    please_select_the_deanery_hospital_you_are_completing_the_form_for, "[\r\n]", " "
  ))) %>%
  left_join(
    volume_clean %>% select(hospital_clean, volume, volume_n),
    by = "hospital_clean"
  )

# Amended to be continuous variable  
#    mutate(hospital_clean = str_squish(str_replace_all(
#    please_select_the_deanery_hospital_you_are_completing_the_form_for, "[\r\n]", " "
#  ))) %>%
#  mutate(
#    volume = case_when(
#      hospital_clean %in% (volume_clean %>% filter(volume == "High") %>% pull(hospital_clean)) ~ "High",
#      hospital_clean %in% (volume_clean %>% filter(volume == "Medium") %>% pull(hospital_clean)) ~ #"Medium",
#      hospital_clean %in% (volume_clean %>% filter(volume == "Low") %>% pull(hospital_clean)) ~ "Low",
#      TRUE ~ NA_character_
#    )
#  )

volume %>% 
  select(volume) %>%
  group_by(volume) %>%
  summarise("Number of Units" = n()) %>%
  gt() %>% 
  gt_theme_espn()

volume	Number of Units
High	1
Medium	8
Low	65

3.18 Location of Injury

injury_location <- snap %>%
  subset(select = c(locality,
                    laterality,
                    where_was_the_site_of_the_fracture_answer_all_that_apply))

predictors <- predictors %>% mutate(
  n_locations = sum(c(
      str_detect(where_was_the_site_of_the_fracture_answer_all_that_apply, "Distal"),
      str_detect(where_was_the_site_of_the_fracture_answer_all_that_apply, "Mid"),
      str_detect(where_was_the_site_of_the_fracture_answer_all_that_apply, "Proximal"),
      str_detect(where_was_the_site_of_the_fracture_answer_all_that_apply, "Crura")
    )),
  proximal = case_when(
    str_detect(where_was_the_site_of_the_fracture_answer_all_that_apply, "Proximal") ~ "Yes",
    TRUE ~ "No"
  ),
  mid = case_when(
    str_detect(where_was_the_site_of_the_fracture_answer_all_that_apply, "Mid") ~ "Yes",
    TRUE ~ "No"
  ),
  distal = case_when(
    str_detect(where_was_the_site_of_the_fracture_answer_all_that_apply, "Distal") ~ "Yes",
    TRUE ~ "No"
  ),
  crura = case_when(
    str_detect(where_was_the_site_of_the_fracture_answer_all_that_apply, "Crura") ~ "Yes",
    TRUE ~ "No"
  ),
  multi_level = case_when(
    n_locations > 1 ~ "Yes",
    TRUE ~ "No"
  ),
  imaging = case_when(
    was_imaging_performed_if_so_what_imaging_select_all_that_apply == "No" ~ "No",
    was_imaging_performed_if_so_what_imaging_select_all_that_apply %in% c("MRI & US",
                                                                          "MRI",
                                                                          "US") ~ "Yes",
    TRUE ~ "No"
    ),
  .keep = "all"
  )  %>% select(-n_locations) %>% cbind()

tab <- predictors %>%
  drop_na(locality, laterality) %>%
  tabyl(locality, laterality) %>%
  adorn_totals("row") %>%
  adorn_totals("col") %>%
  adorn_percentages("row")  

# Get counts separately
tab_counts <- predictors %>%
  drop_na(locality, laterality) %>%
  tabyl(locality, laterality) %>%
  adorn_totals("row") %>%
  adorn_totals("col")

# Combine counts and percentages into n (%)
tab_n_pct <- tab_counts
for(col in 2:ncol(tab_counts)) {
  tab_n_pct[[col]] <- paste0(tab_counts[[col]], " (", 
                             round(tab[[col]]*100, 1), "%)")
}

tab_n_pct %>%
  kable("html", escape = FALSE, 
        col.names = c("Location", "Bilateral", "Unilateral", "Total")) %>%
  kable_styling(full_width = FALSE, 
                bootstrap_options = c("striped", "hover", "condensed")) %>%
  row_spec(0, bold = TRUE, background = "#D3D3D3")

Location	Bilateral	Unilateral	Total
Crura	0 (0%)	9 (100%)	9 (100%)
Distal	1 (4.3%)	22 (95.7%)	23 (100%)
Mid	7 (9.5%)	67 (90.5%)	74 (100%)
Multi-level	3 (42.9%)	4 (57.1%)	7 (100%)
Proximal	14 (12.3%)	100 (87.7%)	114 (100%)
Total	25 (11%)	202 (89%)	227 (100%)

3.19 Final sort of Predictors

predictors_pres_to_th <- predictors %>% 
  subset(select = -c(what_time_did_surgical_repair_take_place_knife_to_skin,
                     what_was_the_estimated_date_and_time_of_injury,
                     when_was_the_patient_first_seen_in_hospital_by_a_healthcare_professional,
                     submission_date,
                     please_select_the_deanery_hospital_you_are_completing_the_form_for,
                     what_suture_was_used_for_the_repair_of_the_tunical_fracture
                     ))

predictors_pres_to_th %>% head() %>% as_tabyl() %>% gt() %>% gt_theme_espn()

submission_id	what_age_was_the_patient_at_the_time_of_the_injury	mechanism_of_injury	if_the_injury_was_sustained_during_sexual_intercourse_what_sexual_position_was_involved	did_the_patient_have_a_history_of_a_prior_penile_fracture	was_pain_present_at_presentation	was_there_immediate_detumescence_loss_of_erection	was_an_audible_pop_reported_by_the_patient	was_penile_bruising_present	did_bruising_extend_over_the_scrotum_lower_abdomen_or_perineum_answer_more_than_one_if_necessary	was_there_a_suspicion_of_urethral_injury_at_presentation	if_there_was_clinical_suspicion_of_a_urethral_injury_why_was_this_select_all_that_apply	was_imaging_performed_if_so_what_imaging_select_all_that_apply	was_the_patient_transferred_between_hospitals_for_specialist_imaging	did_the_reporting_radiologist_have_a_specialist_interest_in_urology	did_the_patient_get_to_theatre_within_24_hours_of_injury	did_the_patient_get_to_theatre_within_24_hours_of_presentation	was_the_patient_transferred_from_another_hospital_for_surgical_repair	were_antibiotics_given_at_induction	was_a_consultant_present_in_theatre	was_a_trainee_present_in_theatre	what_was_the_grade_of_the_primary_surgeon	was_a_specialist_opinion_sought_from_a_urologist_with_a_specialist_interest_in_andrology	was_a_urologist_with_a_specialist_interest_in_andrology_present_in_theatre	where_was_the_site_of_the_fracture_answer_all_that_apply	what_incision_was_made	was_a_concurrent_circumcision_performed	were_intraoperative_urethral_investigations_performed_select_all_that_apply	was_a_urethral_injury_noted_at_surgery	if_there_was_a_urethral_injury_what_suture_was_used_to_repair_the_urethra	was_follow_up_arranged	presentation_to_knife_to_skin	laterality	locality	type_of_suture	any_imaging	hospital_clean	volume	volume_n	proximal	mid	distal	crura	multi_level	imaging
1	53	Sexual intercourse trauma	NA	No	No	NA	NA	Yes	NA	No	No clinical suspicion of urethral injury	US	No	NA	No	No	No	Yes	Yes	Yes	Consultant	NA	NA	Distal shaft right corpora	longitudinal	No	Not performed	No	NA	Yes	527052	Unilateral	Distal	Non-absorbable	Yes	Wales Morriston Hospital, Swansea	Low	1	No	No	Yes	No	No	Yes
2	42	Sexual intercourse trauma	NA	No	Yes	Yes	No	Yes	Scrotum - Lower abdomen	Yes	Blood at urethal meatus	No	No	No	NA	NA	Yes	Yes	Yes	Yes	Specialist Registrar	No	No	Proximal shaft left corpora - Proximal shaft right corpora	Penoscrotal	No	Flexible cystoscopy	Yes	Vicryl Rapide 5-0 Suture	Yes	115	Bilateral	Proximal	Non-absorbable	No	London, South St George's Hospital, London	Medium	6	Yes	No	No	No	No	No
3	44	Sexual intercourse trauma	NA	No	Yes	Yes	NA	No	No bruising present	Yes	Blood at urethal meatus	MRI	No	NA	NA	NA	Yes	Yes	Yes	Yes	Consultant	Yes	Yes	Proximal shaft left corpora - Proximal shaft right corpora	Penoscrotal	No	Not performed	Yes	Vicryl 5-0 Suture	Yes	7620	Bilateral	Proximal	Non-absorbable	Yes	London, South St George's Hospital, London	Medium	6	Yes	No	No	No	No	Yes
4	35	Sexual intercourse trauma	NA	No	Yes	Yes	NA	Yes	No bruising present	Yes	Other	MRI	No	Yes	Yes	Yes	No	Yes	Yes	Yes	Consultant	Yes	Yes	Mid shaft right corpora	longitudinal	NA	Flexible cystoscopy	No	NA	Yes	199	Unilateral	Mid	Non-absorbable	Yes	East Midlands Leicester General Hospital	Low	5	No	Yes	No	No	No	Yes
5	17	NA	NA - Injury not sustained during sexual intercourse	No	Yes	NA	Yes	No	No bruising present	No	No clinical suspicion of urethral injury	No	No	No	No	No	No	Yes	Yes	Yes	Consultant	No	No	NA	Subcoronal with degloving of penis	No	Not performed	No	Other	Yes	2339	NA	NA	Absorbable	No	North West Arrowe Parke Hospital, Wirral	Low	3	No	No	No	No	No	No
6	30	Sexual intercourse trauma	Patient on top (Missionary/face to face)	No	Yes	Yes	Yes	Yes	Scrotum	No	No clinical suspicion of urethral injury	No	No	NA	NA	NA	No	No	Yes	No	Consultant	NA	No	Proximal shaft right corpora	Penoscrotal	No	Not performed	No	Other	Yes	650	Unilateral	Proximal	Unclear	No	Wales Royal Gwent Hospital	Low	2	Yes	No	No	No	No	No

3.20 Plot missing

plot_missing(predictors_pres_to_th)

3.21 Introduce Predictors

plot_bar(predictors_pres_to_th,
         ggtheme = theme_classic(),
         nrow = 1,
         ncol = 1)

## 1 columns ignored with more than 50 categories.
## hospital_clean: 74 categories

4 Explore Outcomes - Descriptive statistics only

4.1 Redefine & Tabulate Outcomes

4.1.1 Redefine

outcomes <- outcomes %>%
  mutate(
    pain = recode_factor(
      at_follow_up_did_the_patient_describe_ongoing_pain,
      "No"  = "No",
      "Yes" = "Yes",
      "No Follow up" = NA_character_,
      "Not from the information available" = NA_character_
    ),
    ed = recode_factor(
      did_the_patient_describe_new_or_worsening_erectile_dysfunction_at_follow_up,
      "No"  = "No",
      "Yes - Manageable with medication" = "Yes",
      "Yes - Still bothersome despite medication" = "Yes",
      "Yes - Unable to quantify from information available" = "Yes",
      "No Follow up" = NA_character_,
      "Not clear from the information available" = NA_character_
    ),
    abnormal_curvature = recode_factor(
      at_follow_up_did_the_patient_describe_penile_curvature,
      "No"  = "No",
      "Yes - Able to have penetrative intercourse" = "Yes",
      "Yes - Not clear if limiting penetrative intercourse" = "Yes",
      "No Follow up" = NA_character_,
      "Not from the information available" = NA_character_
    ),
    waisting_shortening = recode_factor(
      at_follow_up_was_there_evidence_of_penile_wasting_or_shortening,
      "No"  = "No",
      "Yes" = "Yes",
      "No Follow up" = NA_character_,
      "Not clear from the information available" = NA_character_
    ),
    submission_id = submission_id,
    .keep = "unused"
  )

4.1.2 Plot missing

plot_missing(outcomes %>% select(-submission_id))

4.1.3 Tabulate

variable_labels <- c(
  pain = "Ongoing pain",
  ed = "Erectile dysfunction",
  abnormal_curvature = "Abnormal curvature",
  waisting_shortening = "Waisting/shortening"
)

outcomes_summary <- outcomes %>%
  select(-submission_id) %>%
  pivot_longer(cols = c(pain, ed, abnormal_curvature, waisting_shortening),
               names_to = "variable",
               values_to = "value") %>%
  group_by(variable) %>%
  summarise(
    n_total   = sum(!is.na(value)),
    n_missing = sum(is.na(value)),
    n_yes     = sum(value == "Yes", na.rm = TRUE),
    n_no      = sum(value == "No", na.rm = TRUE),
    pct_yes   = round(100 * n_yes / n_total, 1),
    .groups = "drop"
  )

outcomes_summary %>%
  mutate(variable = variable_labels[variable]) %>%
  select(c(variable,
           n_total,
           n_yes,
           pct_yes)) %>%
  gt() %>%
  gt_theme_espn() %>%
  fmt_number(columns = c(pct_yes), decimals = 1) %>%
  cols_label(
    variable   = "Outcome",
    n_total    = "Total (n, non-missing)",
    n_yes      = "Yes (n)",
    pct_yes    = "Yes (%)"
  )

Outcome	Total (n, non-missing)	Yes (n)	Yes (%)
Abnormal curvature	205	40	19.5
Erectile dysfunction	213	39	18.3
Ongoing pain	158	28	17.7
Waisting/shortening	194	17	8.8

4.2 Sort Data prior to Logistic Regression

4.2.1 Pain

predictors_pres_to_th$laterality <- as.factor(predictors_pres_to_th$laterality)
predictors_pres_to_th$locality <- as.factor(predictors_pres_to_th$locality)

predictors_pres_to_th$proximal <- as.factor(predictors_pres_to_th$proximal)
predictors_pres_to_th$distal <- as.factor(predictors_pres_to_th$distal)
predictors_pres_to_th$mid <- as.factor(predictors_pres_to_th$mid)
predictors_pres_to_th$crura <- as.factor(predictors_pres_to_th$crura)
predictors_pres_to_th$multi_level <- as.factor(predictors_pres_to_th$multi_level)
predictors_pres_to_th$imaging <- as.factor(predictors_pres_to_th$imaging)


pain_for_log_reg <- outcomes %>%
  select(pain) %>%
  cbind(predictors_pres_to_th) %>%
  subset(
    select = -c(
      locality,
      if_the_injury_was_sustained_during_sexual_intercourse_what_sexual_position_was_involved,
      if_there_was_clinical_suspicion_of_a_urethral_injury_why_was_this_select_all_that_apply,
      what_was_the_grade_of_the_primary_surgeon,
      were_intraoperative_urethral_investigations_performed_select_all_that_apply,
      was_follow_up_arranged,
      if_there_was_a_urethral_injury_what_suture_was_used_to_repair_the_urethra,
      was_a_specialist_opinion_sought_from_a_urologist_with_a_specialist_interest_in_andrology,
      was_a_concurrent_circumcision_performed,
      was_a_trainee_present_in_theatre,
      did_bruising_extend_over_the_scrotum_lower_abdomen_or_perineum_answer_more_than_one_if_necessary
    )
  ) 

pain_for_log_reg_18 <- pain_for_log_reg %>%
  mutate(pain = pain %>% factor(levels = c("Yes",
                                           "No"))) %>%
  mutate(
    presentation_to_knife_to_skin = presentation_to_knife_to_skin / (60),
    did_the_patient_get_to_theatre_within_18_hours_of_presentation = case_when(
      presentation_to_knife_to_skin >18 ~ "No",
      presentation_to_knife_to_skin <=18 ~ "Yes",
      TRUE ~ NA_character_
    ),
    .keep = "all"
  ) %>% select(-c(did_the_patient_get_to_theatre_within_24_hours_of_presentation,
                  did_the_patient_get_to_theatre_within_24_hours_of_injury))
  
pain_for_log_reg_36 <- pain_for_log_reg %>%
  mutate(pain = pain %>% factor(levels = c("Yes",
                                           "No"))) %>%
  mutate(
    presentation_to_knife_to_skin = presentation_to_knife_to_skin / (60),
    did_the_patient_get_to_theatre_within_36_hours_of_presentation = case_when(
      presentation_to_knife_to_skin >36 ~ "No",
      presentation_to_knife_to_skin <=36 ~ "Yes",
      TRUE ~ NA_character_
    ),
    .keep = "all"
  ) %>% select(-c(did_the_patient_get_to_theatre_within_24_hours_of_presentation,
                  did_the_patient_get_to_theatre_within_24_hours_of_injury)) 

pain_for_log_reg_48 <- pain_for_log_reg %>%
  mutate(pain = pain %>% factor(levels = c("Yes",
                                           "No"))) %>%
  mutate(
    presentation_to_knife_to_skin = presentation_to_knife_to_skin / (60),
    did_the_patient_get_to_theatre_within_48_hours_of_presentation = case_when(
      presentation_to_knife_to_skin >48 ~ "No",
      presentation_to_knife_to_skin <=48 ~ "Yes",
      TRUE ~ NA_character_
    ),
    .keep = "all"
  ) %>% select(-c(did_the_patient_get_to_theatre_within_24_hours_of_presentation,
                  did_the_patient_get_to_theatre_within_24_hours_of_injury)) 

plot_missing(pain_for_log_reg)

glue("Number with complete data: {nrow(pain_for_log_reg %>% drop_na())}")

## Number with complete data: 55

4.2.2 ED

ed_for_log_reg <- outcomes %>%
  select(ed) %>%
  mutate(ed = case_when(ed == "No follow up" ~ NA,
                        TRUE ~ ed) %>% factor(levels = c("Yes",
                                                         "No")),
         .keep = "unused") %>%
  cbind(predictors_pres_to_th) %>%
  subset(select = -c(
    locality,
    if_the_injury_was_sustained_during_sexual_intercourse_what_sexual_position_was_involved,
    if_there_was_clinical_suspicion_of_a_urethral_injury_why_was_this_select_all_that_apply,
    what_was_the_grade_of_the_primary_surgeon,
    were_intraoperative_urethral_investigations_performed_select_all_that_apply,
    was_follow_up_arranged,
    if_there_was_a_urethral_injury_what_suture_was_used_to_repair_the_urethra,
    was_a_specialist_opinion_sought_from_a_urologist_with_a_specialist_interest_in_andrology,
    was_a_concurrent_circumcision_performed,
    was_a_trainee_present_in_theatre,
    did_bruising_extend_over_the_scrotum_lower_abdomen_or_perineum_answer_more_than_one_if_necessary
  )) 

ed_for_log_reg_18 <- ed_for_log_reg %>%
  mutate(
    presentation_to_knife_to_skin = presentation_to_knife_to_skin / (60),
    did_the_patient_get_to_theatre_within_18_hours_of_presentation = case_when(
      presentation_to_knife_to_skin >18 ~ "No",
      presentation_to_knife_to_skin <=18 ~ "Yes",
      TRUE ~ NA_character_
    ),
    .keep = "all"
  ) %>% select(-c(did_the_patient_get_to_theatre_within_24_hours_of_presentation,
                  did_the_patient_get_to_theatre_within_24_hours_of_injury))
  
ed_for_log_reg_36 <- ed_for_log_reg %>%
  mutate(
    presentation_to_knife_to_skin = presentation_to_knife_to_skin / (60),
    did_the_patient_get_to_theatre_within_36_hours_of_presentation = case_when(
      presentation_to_knife_to_skin >36 ~ "No",
      presentation_to_knife_to_skin <=36 ~ "Yes",
      TRUE ~ NA_character_
    ),
    .keep = "all"
  ) %>% select(-c(did_the_patient_get_to_theatre_within_24_hours_of_presentation,
                  did_the_patient_get_to_theatre_within_24_hours_of_injury)) 

ed_for_log_reg_48 <- ed_for_log_reg %>%
  mutate(
    presentation_to_knife_to_skin = presentation_to_knife_to_skin / (60),
    did_the_patient_get_to_theatre_within_48_hours_of_presentation = case_when(
      presentation_to_knife_to_skin >48 ~ "No",
      presentation_to_knife_to_skin <=48 ~ "Yes",
      TRUE ~ NA_character_
    ),
    .keep = "all"
  ) %>% select(-c(did_the_patient_get_to_theatre_within_24_hours_of_presentation,
                  did_the_patient_get_to_theatre_within_24_hours_of_injury)) 

plot_missing(ed_for_log_reg)

glue("Number with complete data: {nrow(ed_for_log_reg %>% drop_na())}")

## Number with complete data: 52

4.2.3 Waisting/Shortening

waisting_shortening_for_log_reg <- outcomes %>%
  select(waisting_shortening) %>%
  mutate(waisting_shortening = case_when(waisting_shortening == "No follow up" ~ NA,
                        TRUE ~ waisting_shortening) %>% factor(levels = c("Yes",
                                                         "No")),
         .keep = "unused") %>%
  cbind(predictors_pres_to_th) %>%
  subset(select = -c(
    if_the_injury_was_sustained_during_sexual_intercourse_what_sexual_position_was_involved,
    if_there_was_clinical_suspicion_of_a_urethral_injury_why_was_this_select_all_that_apply,
    what_was_the_grade_of_the_primary_surgeon,
    were_intraoperative_urethral_investigations_performed_select_all_that_apply,
    was_follow_up_arranged,
    if_there_was_a_urethral_injury_what_suture_was_used_to_repair_the_urethra,
    was_a_specialist_opinion_sought_from_a_urologist_with_a_specialist_interest_in_andrology,
    was_a_concurrent_circumcision_performed,
    was_a_trainee_present_in_theatre,
    did_bruising_extend_over_the_scrotum_lower_abdomen_or_perineum_answer_more_than_one_if_necessary
  )) 

waisting_shortening_for_log_reg_18 <- waisting_shortening_for_log_reg %>%
  mutate(
    presentation_to_knife_to_skin = presentation_to_knife_to_skin / (60),
    did_the_patient_get_to_theatre_within_18_hours_of_presentation = case_when(
      presentation_to_knife_to_skin >18 ~ "No",
      presentation_to_knife_to_skin <=18 ~ "Yes",
      TRUE ~ NA_character_
    ),
    .keep = "all"
  ) %>% select(-c(did_the_patient_get_to_theatre_within_24_hours_of_presentation,
                  did_the_patient_get_to_theatre_within_24_hours_of_injury))
  
waisting_shortening_for_log_reg_36 <- waisting_shortening_for_log_reg %>%
  mutate(
    presentation_to_knife_to_skin = presentation_to_knife_to_skin / (60),
    did_the_patient_get_to_theatre_within_36_hours_of_presentation = case_when(
      presentation_to_knife_to_skin >36 ~ "No",
      presentation_to_knife_to_skin <=36 ~ "Yes",
      TRUE ~ NA_character_
    ),
    .keep = "all"
  ) %>% select(-c(did_the_patient_get_to_theatre_within_24_hours_of_presentation,
                  did_the_patient_get_to_theatre_within_24_hours_of_injury)) 

waisting_shortening_for_log_reg_48 <- waisting_shortening_for_log_reg %>%
  mutate(
    presentation_to_knife_to_skin = presentation_to_knife_to_skin / (60),
    did_the_patient_get_to_theatre_within_48_hours_of_presentation = case_when(
      presentation_to_knife_to_skin >48 ~ "No",
      presentation_to_knife_to_skin <=48 ~ "Yes",
      TRUE ~ NA_character_
    ),
    .keep = "all"
  ) %>% select(-c(did_the_patient_get_to_theatre_within_24_hours_of_presentation,
                  did_the_patient_get_to_theatre_within_24_hours_of_injury)) 

plot_missing(waisting_shortening_for_log_reg)

glue("Number with complete data: {nrow(waisting_shortening_for_log_reg %>% drop_na())}")

## Number with complete data: 51

4.2.4 Curvature

abnormal_curvature_for_log_reg <- outcomes %>%
  select(abnormal_curvature) %>%
  mutate(abnormal_curvature = case_when(abnormal_curvature == "No follow up" ~ NA,
                        TRUE ~ abnormal_curvature) %>% factor(levels = c("Yes",
                                                         "No")),
         .keep = "unused") %>%
  cbind(predictors_pres_to_th) %>%
  subset(select = -c(
    locality,
    if_the_injury_was_sustained_during_sexual_intercourse_what_sexual_position_was_involved,
    if_there_was_clinical_suspicion_of_a_urethral_injury_why_was_this_select_all_that_apply,
    what_was_the_grade_of_the_primary_surgeon,
    were_intraoperative_urethral_investigations_performed_select_all_that_apply,
    was_follow_up_arranged,
    if_there_was_a_urethral_injury_what_suture_was_used_to_repair_the_urethra,
    was_a_specialist_opinion_sought_from_a_urologist_with_a_specialist_interest_in_andrology,
    was_a_concurrent_circumcision_performed,
    was_a_trainee_present_in_theatre,
    did_bruising_extend_over_the_scrotum_lower_abdomen_or_perineum_answer_more_than_one_if_necessary
  )) 

abnormal_curvature_for_log_reg_18 <- abnormal_curvature_for_log_reg %>%
  mutate(
    presentation_to_knife_to_skin = presentation_to_knife_to_skin / (60),
    did_the_patient_get_to_theatre_within_18_hours_of_presentation = case_when(
      presentation_to_knife_to_skin >18 ~ "No",
      presentation_to_knife_to_skin <=18 ~ "Yes",
      TRUE ~ NA_character_
    ),
    .keep = "all"
  ) %>% select(-c(did_the_patient_get_to_theatre_within_24_hours_of_presentation,
                  did_the_patient_get_to_theatre_within_24_hours_of_injury))
  
abnormal_curvature_for_log_reg_36 <- abnormal_curvature_for_log_reg %>%
  mutate(
    presentation_to_knife_to_skin = presentation_to_knife_to_skin / (60),
    did_the_patient_get_to_theatre_within_36_hours_of_presentation = case_when(
      presentation_to_knife_to_skin >36 ~ "No",
      presentation_to_knife_to_skin <=36 ~ "Yes",
      TRUE ~ NA_character_
    ),
    .keep = "all"
  ) %>% select(-c(did_the_patient_get_to_theatre_within_24_hours_of_presentation,
                  did_the_patient_get_to_theatre_within_24_hours_of_injury)) 

abnormal_curvature_for_log_reg_48 <- abnormal_curvature_for_log_reg %>%
  mutate(
    presentation_to_knife_to_skin = presentation_to_knife_to_skin / (60),
    did_the_patient_get_to_theatre_within_48_hours_of_presentation = case_when(
      presentation_to_knife_to_skin >48 ~ "No",
      presentation_to_knife_to_skin <=48 ~ "Yes",
      TRUE ~ NA_character_
    ),
    .keep = "all"
  ) %>% select(-c(did_the_patient_get_to_theatre_within_24_hours_of_presentation,
                  did_the_patient_get_to_theatre_within_24_hours_of_injury)) 

plot_missing(abnormal_curvature_for_log_reg)

glue("Number with complete data: {nrow(abnormal_curvature_for_log_reg %>% drop_na())}")

## Number with complete data: 50

4.3 Explore Time to Theatre Effect on Outcomes

4.3.1 Pain

4.3.1.1 Histogram

time_from_injury_to_th1 <- time_from_injury_to_th %>% 
  subset(select = c(
    submission_id,
    time_from_injury_to_th
  ))

knife_to_skin <- predictors %>% 
  select(submission_id,
         presentation_to_knife_to_skin) %>%
  left_join(outcomes, by = c("submission_id" = "submission_id")) %>% 
  left_join(time_from_injury_to_th1, by = c("submission_id" = "submission_id"))

p1 <- knife_to_skin %>% 
  drop_na(presentation_to_knife_to_skin, pain) %>%
  ggplot(aes(x = presentation_to_knife_to_skin / 60, fill = pain)) +
  geom_histogram(bins = 50) +
  labs(
    title = "Histogram of Time from Presentation to Theatre",
    x = "Hours",
    y = "Count"
  ) +
  theme_minimal() + xlim(0,100)

p2 <- ggplot(knife_to_skin %>% drop_na(time_from_injury_to_th,
                                       pain)
             , aes(x = time_from_injury_to_th / 60, fill = pain)) +
  geom_histogram(bins = 50) +
  labs(
    title = "Histogram of Time from Injury to Theatre",
    x = "Hours",
    y = "Count"
  ) +
  theme_minimal()

cowplot::plot_grid(p1, p2)

## Warning: Removed 20 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 4 rows containing missing values or values outside the scale range
## (`geom_bar()`).

4.3.1.2 Density plots

p1 <- ggplot(knife_to_skin %>% drop_na(presentation_to_knife_to_skin,
                                       pain),
             aes(x = presentation_to_knife_to_skin / 60, fill = pain)) +
  geom_density() +
  labs(
    title = "Time from Presentation to Theatre",
    x = "Hours",
    y = "Count"
  ) + xlim(0,100) +
  theme_minimal()

p2 <- ggplot(knife_to_skin %>% drop_na(time_from_injury_to_th,
                                       pain)
             , aes(x = time_from_injury_to_th / 60, fill = pain)) +
  geom_density() +
  labs(
    title = "Time from Injury to Theatre",
    x = "Hours",
    y = "Count"
  ) + xlim(0,100) +
  theme_minimal()

cowplot::plot_grid(p1, p2)

## Warning: Removed 20 rows containing non-finite outside the scale range
## (`stat_density()`).

## Warning: Removed 7 rows containing non-finite outside the scale range
## (`stat_density()`).

4.3.1.3 Bins

time_from_injury_to_th1 <- time_from_injury_to_th %>% 
  subset(select = c(
    submission_id,
    time_from_injury_to_th_bins,
    presentation_to_knife_to_skin_bins
  ))

knife_to_skin1 <- time_from_injury_to_th1 %>% 
  left_join(outcomes, by = c("submission_id" = "submission_id")) 

p1 <- knife_to_skin1 %>% 
  drop_na(presentation_to_knife_to_skin_bins) %>%
  group_by(presentation_to_knife_to_skin_bins) %>%
  summarise(
    n_yes = sum(pain == "Yes", na.rm = TRUE),
    n_total = n(),
    percent_yes = 100 * n_yes / n_total,
    .groups = "drop"
  ) %>%
  ggplot(aes(x = presentation_to_knife_to_skin_bins, y = percent_yes)) +
  geom_col(fill = "#D7191C") +
  geom_text(aes(label = paste0(round(percent_yes, 1), "%")),
            vjust = -0.5, size = 3.5) +
  labs(
    title = "Percentage of Patients Reporting Pain",
    x = "Time to Theatre from Presentation (Hours, binned)",
    y = "Percentage with Pain"
  ) +
  theme_minimal()

p2 <- knife_to_skin1 %>% 
  drop_na(presentation_to_knife_to_skin_bins) %>%
  group_by(presentation_to_knife_to_skin_bins) %>%
  summarise(
    n_yes = sum(pain == "Yes", na.rm = TRUE),
    n_total = n(),
    percent_yes = 100 * n_yes / n_total,
    .groups = "drop"
  ) %>%
  ggplot(aes(x = presentation_to_knife_to_skin_bins, y = n_yes)) +
  geom_col(fill = "#D7191C") +
  geom_text(aes(label = paste0("n=",round(n_yes, 1))),
            vjust = -0.5, size = 3.5) +
  labs(
    title = "Number of Patients Reporting Pain",
    x = "Time to Theatre from Presentation (Hours, binned)",
    y = "Number with Pain"
  ) +
  theme_minimal()


p3 <- knife_to_skin1 %>% 
  drop_na(time_from_injury_to_th_bins) %>%
  group_by(time_from_injury_to_th_bins) %>%
  summarise(
    n_yes = sum(pain == "Yes", na.rm = TRUE),
    n_total = n(),
    percent_yes = 100 * n_yes / n_total,
    .groups = "drop"
  ) %>%
  ggplot(aes(x = time_from_injury_to_th_bins, y = percent_yes)) +
  geom_col(fill = "#D7191C") +
  geom_text(aes(label = paste0(round(percent_yes, 1), "%")),
            vjust = -0.5, size = 3.5) +
  labs(
    title = "Percentage of Patients Reporting Pain",
    x = "Time to Theatre from Injury (Hours, binned)",
    y = "Percentage with Pain"
  ) +
  theme_minimal()

p4 <- knife_to_skin1 %>% 
  drop_na(time_from_injury_to_th_bins) %>%
  group_by(time_from_injury_to_th_bins) %>%
  summarise(
    n_yes = sum(pain == "Yes", na.rm = TRUE),
    n_total = n(),
    percent_yes = 100 * n_yes / n_total,
    .groups = "drop"
  ) %>%
  ggplot(aes(x = time_from_injury_to_th_bins, y = n_yes)) +
  geom_col(fill = "#D7191C") +
  geom_text(aes(label = paste0("n=", round(n_yes, 1))),
            vjust = -0.5, size = 3.5) +
  labs(
    title = "Number of Patients Reporting Pain",
    x = "Time to Theatre from Injury (Hours, binned)",
    y = "Number with Pain"
  ) +
  theme_minimal()

cowplot::plot_grid(p1, p2)

cowplot::plot_grid(p3, p4)

4.3.2 ED

4.3.2.1 Histogram

p1 <- ggplot(knife_to_skin %>% drop_na(presentation_to_knife_to_skin,
                                       ed) %>%
               subset(ed != "No follow up"),
             aes(x = presentation_to_knife_to_skin/60, fill = ed)) +
  geom_histogram(bins = 50) +
  labs(
    title = "Time from Presentation to Theatre",
    x = "Hours",
    y = "Count"
  ) +
  theme_minimal() + xlim(0,100)

p2 <- ggplot(knife_to_skin %>% drop_na(time_from_injury_to_th,
                                       ed) %>%
               subset(ed != "No follow up")
             , aes(x = time_from_injury_to_th/60, fill = ed)) +
  geom_histogram(bins = 50) +
  labs(
    title = "Time from Injury to Theatre",
    x = "Hours",
    y = "Count"
  ) +
  theme_minimal()  + xlim(0,100)

cowplot::plot_grid(p1, p2)

## Warning: Removed 19 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 4 rows containing missing values or values outside the scale range
## (`geom_bar()`).

## Warning: Removed 7 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 4 rows containing missing values or values outside the scale range
## (`geom_bar()`).

4.3.2.2 Density plot

p1 <- ggplot(knife_to_skin %>% drop_na(presentation_to_knife_to_skin,
                                       ed) %>% subset(
                                         ed != "No follow up"
                                       ),
             aes(x = presentation_to_knife_to_skin/60, fill = ed)) +
  geom_density() +
  labs(
    title = "Time from Presentation to Theatre",
    x = "Hours",
    y = "Count"
  ) +
  theme_minimal()  + xlim(0,100)

p2 <- ggplot(knife_to_skin %>% drop_na(time_from_injury_to_th,
                                       ed) %>% subset(
                                         ed != "No follow up"
                                       )
             , aes(x = time_from_injury_to_th/60, fill = ed)) +
  geom_density() +
  labs(
    title = "Time from Injury to Theatre",
    x = "Hours",
    y = "Count"
  ) +
  theme_minimal()  + xlim(0,100)

cowplot::plot_grid(p1, p2)

## Warning: Removed 19 rows containing non-finite outside the scale range
## (`stat_density()`).

## Warning: Removed 7 rows containing non-finite outside the scale range
## (`stat_density()`).

4.3.2.3 Bins

time_from_injury_to_th1 <- time_from_injury_to_th %>% 
  subset(select = c(
    submission_id,
    time_from_injury_to_th_bins,
    presentation_to_knife_to_skin_bins
  ))

knife_to_skin1 <- time_from_injury_to_th1 %>% 
  left_join(outcomes, by = c("submission_id" = "submission_id")) 

p1 <- knife_to_skin1 %>% 
  drop_na(presentation_to_knife_to_skin_bins,
          ed) %>%
  group_by(presentation_to_knife_to_skin_bins) %>%
  summarise(
    n_yes = sum(ed == "Yes", na.rm = TRUE),
    n_total = n(),
    percent_yes = 100 * n_yes / n_total,
    .groups = "drop"
  ) %>%
  ggplot(aes(x = presentation_to_knife_to_skin_bins, y = percent_yes)) +
  geom_col(fill = "#D7191C") +
  geom_text(aes(label = paste0(round(percent_yes, 1), "%")),
            vjust = -0.5, size = 3.5) +
  labs(
    title = "Percentage of Patients Reporting ED",
    x = "Time to Theatre from Presentation (Hours, binned)",
    y = "Percentage with ED"
  ) +
  theme_minimal()

p2 <- knife_to_skin1 %>% 
  drop_na(presentation_to_knife_to_skin_bins,
          ed) %>%
  group_by(presentation_to_knife_to_skin_bins) %>%
  summarise(
    n_yes = sum(ed == "Yes", na.rm = TRUE),
    n_total = n(),
    percent_yes = 100 * n_yes / n_total,
    .groups = "drop"
  ) %>%
  ggplot(aes(x = presentation_to_knife_to_skin_bins, y = n_yes)) +
  geom_col(fill = "#D7191C") +
   geom_text(aes(label = paste0("n=",round(n_yes, 1))),
            vjust = -0.5, size = 3.5) +
  labs(
    title = "Number of Patients Reporting ED",
    x = "Time to Theatre from Presentation (Hours, binned)",
    y = "Number with ED"
  ) +
  theme_minimal()


p3 <- knife_to_skin1 %>% 
  drop_na(time_from_injury_to_th_bins,
          ed) %>%
  group_by(time_from_injury_to_th_bins) %>%
  summarise(
    n_yes = sum(ed == "Yes", na.rm = TRUE),
    n_total = n(),
    percent_yes = 100 * n_yes / n_total,
    .groups = "drop"
  ) %>%
  ggplot(aes(x = time_from_injury_to_th_bins, y = percent_yes)) +
  geom_col(fill = "#D7191C") +
  geom_text(aes(label = paste0(round(percent_yes, 1), "%")),
            vjust = -0.5, size = 3.5) +
  labs(
    title = "Percentage of Patients Reporting ED",
    x = "Time to Theatre from Injury (Hours, binned)",
    y = "Percentage with ED"
  ) +
  theme_minimal()

p4 <- knife_to_skin1 %>% 
  drop_na(time_from_injury_to_th_bins,
          ed) %>%
  group_by(time_from_injury_to_th_bins) %>%
  summarise(
    n_yes = sum(ed == "Yes", na.rm = TRUE),
    n_total = n(),
    percent_yes = 100 * n_yes / n_total,
    .groups = "drop"
  ) %>%
  ggplot(aes(x = time_from_injury_to_th_bins, y = n_yes)) +
  geom_col(fill = "#D7191C") +
  geom_text(aes(label = paste0("n=",round(n_yes, 1))),
            vjust = -0.5, size = 3.5) +
  labs(
    title = "Number of Patients Reporting ED",
    x = "Time to Theatre from Injury (Hours, binned)",
    y = "Number with ED"
  ) +
  theme_minimal()

cowplot::plot_grid(p1, p2)

cowplot::plot_grid(p3, p4)

4.3.3 Waisting/Shortening

4.3.3.1 Histogram

p1 <- ggplot(knife_to_skin %>% drop_na(presentation_to_knife_to_skin,
                                       waisting_shortening) %>%
               subset(waisting_shortening != "No follow up"),
             aes(x = presentation_to_knife_to_skin/60, fill = waisting_shortening)) +
  geom_histogram(bins = 50) +
  labs(
    title = "Time from Presentation to Theatre",
    x = "Hours",
    y = "Count"
  ) +
  theme_minimal() +xlim(0,100)

p2 <- ggplot(knife_to_skin %>% drop_na(time_from_injury_to_th,
                                       waisting_shortening) %>%
               subset(waisting_shortening != "No follow up")
             , aes(x = time_from_injury_to_th/60, fill = waisting_shortening)) +
  geom_histogram(bins = 50) +
  labs(
    title = "Time from Injury to Theatre",
    x = "Hours",
    y = "Count"
  ) +
  theme_minimal()

cowplot::plot_grid(p1, p2)

## Warning: Removed 15 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 4 rows containing missing values or values outside the scale range
## (`geom_bar()`).

4.3.3.2 Density plot

p1 <- ggplot(knife_to_skin %>% drop_na(presentation_to_knife_to_skin,
                                       waisting_shortening) %>% subset(
                                         waisting_shortening != "No follow up"
                                       ),
             aes(x = presentation_to_knife_to_skin/60, fill = waisting_shortening)) +
  geom_density() +
  labs(
    title = "Time from Presentation to Theatre",
    x = "Hours",
    y = "Count"
  ) +
  theme_minimal() + xlim(0,100)

p2 <- ggplot(knife_to_skin %>% drop_na(time_from_injury_to_th,
                                       waisting_shortening) %>% subset(
                                         waisting_shortening != "No follow up"
                                       )
             , aes(x = time_from_injury_to_th/60, fill = waisting_shortening)) +
  geom_density() +
  labs(
    title = "Time from Injury to Theatre",
    x = "Hours",
    y = "Count"
  ) +
  theme_minimal()

cowplot::plot_grid(p1, p2)

## Warning: Removed 15 rows containing non-finite outside the scale range
## (`stat_density()`).

4.3.4 Curvature

4.3.4.1 Histogram

p1 <- ggplot(knife_to_skin %>% drop_na(presentation_to_knife_to_skin,
                                       abnormal_curvature) %>%
               subset(abnormal_curvature != "No follow up"),
             aes(x = presentation_to_knife_to_skin/60, fill = abnormal_curvature)) +
  geom_histogram(bins = 50) +
  labs(
    title = "Time from Presentation to Theatre",
    x = "Hours",
    y = "Count"
  ) +
  theme_minimal() + xlim(0,100)

p2 <- ggplot(knife_to_skin %>% drop_na(time_from_injury_to_th,
                                       abnormal_curvature) %>%
               subset(abnormal_curvature != "No follow up")
             , aes(x = time_from_injury_to_th/60, fill = abnormal_curvature)) +
  geom_histogram(bins = 50) +
  labs(
    title = "Time from Injury to Theatre",
    x = "Hours",
    y = "Count"
  ) +
  theme_minimal()

cowplot::plot_grid(p1, p2)

## Warning: Removed 19 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 4 rows containing missing values or values outside the scale range
## (`geom_bar()`).

4.3.4.2 Density plot

p1 <- ggplot(knife_to_skin %>% drop_na(presentation_to_knife_to_skin,
                                       abnormal_curvature) %>% subset(
                                         abnormal_curvature != "No follow up"
                                       ),
             aes(x = presentation_to_knife_to_skin/60, fill = abnormal_curvature)) +
  geom_density() +
  labs(
    title = "Time from Presentation to Theatre",
    x = "Hours",
    y = "Count"
  ) +
  theme_minimal()  + xlim(0,100)

p2 <- ggplot(knife_to_skin %>% drop_na(time_from_injury_to_th,
                                       abnormal_curvature) %>% subset(
                                         abnormal_curvature != "No follow up"
                                       )
             , aes(x = time_from_injury_to_th/60, fill = abnormal_curvature)) +
  geom_density() +
  labs(
    title = "Time from Injury to Theatre",
    x = "Hours",
    y = "Count"
  ) +
  theme_minimal()  + xlim(0,100)

cowplot::plot_grid(p1, p2)

## Warning: Removed 19 rows containing non-finite outside the scale range
## (`stat_density()`).

## Warning: Removed 6 rows containing non-finite outside the scale range
## (`stat_density()`).

5 Logistic Regression - Time from Presentation to Theatre with 24 hour cut-off

5.1 Functions

5.1.1 Post-Hoc Power analysis Function

# Function to calculate power for two proportions and output as a gt table
calculate_power_table <- function(p1, p2, n1, n2, alpha = 0.05, title = "Power Calculation Results") {
  # Complementary probabilities
  q1 <- 1 - p1
  q2 <- 1 - p2
  
  # Combined proportions
  k <- n1 / n2
  p_bar <- (p1 + k * p2) / (1 + k)
  q_bar <- 1 - p_bar
  
  # Standard error terms
  se_diff <- sqrt((p1 * q1 / n1) + (p2 * q2 / n2))
  se_null <- sqrt(p_bar * q_bar * (1 / n1 + 1 / n2))
  
  # Z-score for significance level (two-sided)
  z_alpha <- qnorm(1 - alpha / 2)
  
  # Z-score under the alternative hypothesis
  z_alt <- (p2 - p1) / se_diff - z_alpha * (se_null / se_diff)
  
  # Calculate power
  power <- pnorm(z_alt)
  
  # Prepare results for table
  results <- data.frame(
    Parameter = c("Group 1 Proportion (p1)", "Group 2 Proportion (p2)",
                  "Group 1 Size (n1)", "Group 2 Size (n2)",
                  "Significance Level (alpha)", "Calculated Power (%)"),
    Value = c(p1, p2, n1, n2, alpha, round(power * 100, 1))
  )
  
  # Create gt table
  power_table <- gt(results) %>%
    tab_header(
      title = title,
      subtitle = "Based on Two Proportion Test"
    ) %>%
    cols_label(
      Parameter = "Parameter",
      Value = "Value"
    ) %>%
    fmt_number(
      columns = "Value",
      decimals = 2,
      rows = 1:4
    )
  
  return(power_table)
}

5.2 Pain

5.2.1 Post-hoc Power calculation

n1 <- nrow(pain_for_log_reg %>% subset(did_the_patient_get_to_theatre_within_24_hours_of_presentation == "No"))
p1 <- nrow(pain_for_log_reg %>% subset(did_the_patient_get_to_theatre_within_24_hours_of_presentation == "No" & pain == "Yes")) / n1

n2 <- nrow(pain_for_log_reg %>% subset(did_the_patient_get_to_theatre_within_24_hours_of_presentation == "Yes"))
alpha <- 0.05

p2 <- nrow(pain_for_log_reg %>% subset(did_the_patient_get_to_theatre_within_24_hours_of_presentation == "Yes" & pain == "Yes")) / n2

power <- calculate_power_table(p1, p2, n1, n2, alpha, title = "Power Calculation Results for Post-Operative Pain")
power

Power Calculation Results for Post-Operative Pain
Based on Two Proportion Test
Parameter	Value
Group 1 Proportion (p1)	0.12
Group 2 Proportion (p2)	0.11
Group 1 Size (n1)	76.00
Group 2 Size (n2)	138.00
Significance Level (alpha)	0.05
Calculated Power (%)	1.50

5.2.2 Logistic Regression (Firth Correction)

knife_to_skin1 <- knife_to_skin %>% 
  subset(select = c(
    submission_id,
    time_from_injury_to_th
  ))

pain_for_log_reg <- pain_for_log_reg %>% 
  left_join(knife_to_skin1,
            by = c("submission_id" = "submission_id")) %>%
  mutate(pain = pain %>% factor(levels = c("Yes",
                                           "No")))

pain_for_log_reg2 <- pain_for_log_reg %>%
  mutate(presentation_to_knife_to_skin = ifelse(
    presentation_to_knife_to_skin < 0 | presentation_to_knife_to_skin > 7920,
    NA,
    presentation_to_knife_to_skin / (60*6)
  ),
    pain = case_when(pain == "Yes" ~ 1,
                     pain == "No" ~ 0,
                     TRUE ~ NA_real_),
  .keep = "unused"
  ) %>%
  select(
    pain,
    imaging,
    type_of_suture,
    volume_n,
    what_age_was_the_patient_at_the_time_of_the_injury,
    did_the_patient_get_to_theatre_within_24_hours_of_presentation,
    what_age_was_the_patient_at_the_time_of_the_injury,
    was_the_patient_transferred_from_another_hospital_for_surgical_repair,
    was_a_urethral_injury_noted_at_surgery,
    laterality,
    distal,
    proximal,
    mid,
    crura,
    multi_level,
    mechanism_of_injury
  ) %>% drop_na()

pain_lr <- logistf(
  pain ~ did_the_patient_get_to_theatre_within_24_hours_of_presentation + imaging + volume_n + type_of_suture +
    what_age_was_the_patient_at_the_time_of_the_injury + was_a_urethral_injury_noted_at_surgery +
    was_the_patient_transferred_from_another_hospital_for_surgical_repair + distal + proximal + mid +
    crura + multi_level + laterality + mechanism_of_injury,
  data = pain_for_log_reg2
)

pain_lr %>%
  tbl_regression(
    exponentiate = TRUE,
    label = list(
      did_the_patient_get_to_theatre_within_24_hours_of_presentation = "Theatre within 24 hours of Presentation",
      was_a_urethral_injury_noted_at_surgery = "Urethral Injury",
      was_the_patient_transferred_from_another_hospital_for_surgical_repair = "Transfer for Surgical Repair",
      laterality = "Side of Injury",
      distal = "Distal Corpus Location",
    proximal = "Proximal Corpus Location",
    mid = "Mid Corpus Location",
    crura = "Crural location",
    multi_level = "Multi-level Injury",
      mechanism_of_injury = "Mechanism of Injury",
      what_age_was_the_patient_at_the_time_of_the_injury =
        "Age at Time of Injury",
    imaging = "Imaging",
    volume_n = "Centre Volume",
    type_of_suture = "Type of Suture"
    )
  ) %>%
  bold_labels() %>%
  as_gt() %>%  
  tab_source_note(
    md(glue(
      "Number of patients in analysis: {nrow(pain_for_log_reg2)}"
    ))
  )

## ! `broom::tidy()` failed to tidy the model.

## ✔ `tidy_parameters()` used instead.

## ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.

Characteristic	exp(Beta)	95% CI	p-value
Theatre within 24 hours of Presentation
No	—	—
Yes	1.70	0.51, 6.42	0.4
Imaging
No	—	—
Yes	1.17	0.33, 4.56	0.8
Centre Volume	0.93	0.87, 0.98	0.002
Type of Suture
Absorbable	—	—
Non-absorbable	0.46	0.10, 2.84	0.4
Unclear	1.12	0.13, 10.6	>0.9
Age at Time of Injury	1.04	0.99, 1.10	0.12
Urethral Injury
No	—	—
Yes	1.69	0.40, 6.37	0.5
Transfer for Surgical Repair
No	—	—
Yes	6.16	1.70, 27.2	0.005
Distal Corpus Location
No	—	—
Yes	0.00	0.00, 8.44	0.2
Proximal Corpus Location
No	—	—
Yes	0.00	0.00, 5.54	0.13
Mid Corpus Location
No	—	—
Yes	0.00	0.00, 2.66	0.086
Crural location
No	—	—
Yes	0.04	0.00, 27.6	0.3
Multi-level Injury
No	—	—
Yes	4,446	0.41, 904,601,017	0.077
Side of Injury
Bilateral	—	—
Unilateral	0.68	0.13, 3.98	0.6
Mechanism of Injury
Sexual intercourse trauma	—	—
Other	0.79	0.01, 13.3	0.9
Non-sexual blunt force trauma	4.35	0.85, 21.3	0.076
Masturbation injury	5.25	0.64, 58.0	0.12
Abbreviation: CI = Confidence Interval
Number of patients in analysis: 116

# Add back in to get 24hr cut-off: presentation_to_knife_to_skin,
# did_the_patient_get_to_theatre_within_24_hours_of_presentation
# time_from_injury_to_th

5.2.3 ROC

pain_lr_predict <-
  predict(pain_lr, pain_for_log_reg, type = "response")

pain_outcome <- as.factor(ifelse(pain_for_log_reg$pain == "Yes",
                                  1,
                                  0))

pain_lr_roc <-
  pROC::roc(pain_outcome, pain_lr_predict)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

auc(pain_lr_roc)

## Area under the curve: 0.8406

ci.auc(pain_lr_roc)

## 95% CI: 0.7538-0.9274 (DeLong)

pain_lr_roc_data <- cbind(sensitivities = pain_lr_roc$sensitivities,
                                specificities = pain_lr_roc$specificities) %>% as_tibble()

ggplot(as_tibble(pain_lr_roc_data), aes(x = 1 - specificities, y = sensitivities)) +
  geom_line() +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed")

5.2.4 Calibration plot

pain_lr_predict2 <- as_tibble(pain_lr_predict) %>% cbind(pain_for_log_reg) %>% subset(select = c(
  value,
  pain,
  did_the_patient_get_to_theatre_within_24_hours_of_presentation
)) %>% mutate(
  pred = value,
  .keep = "unused"
)

pain_lr_predict2$pain <- ifelse(pain_lr_predict2$pain == "Yes", 1, 0)

pain_lr_predict2$did_the_patient_get_to_theatre_within_24_hours_of_presentation <- as.integer(pain_lr_predict2$did_the_patient_get_to_theatre_within_24_hours_of_presentation)

predtools::calibration_plot(
  data = pain_lr_predict2,
  obs = "pain",
  pred = "pred",
  y_lim = c(0, 1),
  x_lim = c(0, 1),
  title = "Calibration plot for Pain Model"
)

## $calibration_plot

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

5.2.5 Brier Score

brier <- function(model, data, outcome) {
  probs <- predict(model, type = "response")
  obs <- data[[outcome]]
  mean((obs - probs)^2, na.rm = TRUE)
}

brier_score <- brier(pain_lr, pain_for_log_reg2, "pain") 

cbind(
  "Score" = "Brier Score",
  "Brier Score" = round(brier_score, digits = 3)
) %>% as_tibble() %>% gt() %>% gt_theme_espn()

Score	Brier Score
Brier Score	0.11

5.2.6 Nagelkerke’s R squared

pseudoR2_logistf <- function(model, data) {
  # Extract log-likelihoods
  ll_full <- model$loglik["full"]
  ll_null <- model$loglik["null"]
  
  # Sample size
  n <- nrow(data)
  
  # Cox & Snell R²
  R2_CS <- 1 - exp((2 / n) * (ll_null - ll_full))
  
  # Nagelkerke R²
  R2_N <- R2_CS / (1 - exp((2 / n) * ll_null))

  R2_N <- R2_N %>% as_tibble()
  
  x <- cbind(
  "Score" = "Nagelkerke's R^2",
  "Nagelkerke's R^2" = round(R2_N$value, digits = 3)
) %>% as_tibble() %>% gt() %>% gt_theme_espn()
  
  # Return as a named list
  return(
    x
  )
}

pseudoR2_logistf(pain_lr, pain_for_log_reg2)

Score	Nagelkerke's R^2
Nagelkerke's R^2	0.365

5.3 ED

5.3.1 Post-hoc Power calculation

n1 <- nrow(ed_for_log_reg %>% subset(did_the_patient_get_to_theatre_within_24_hours_of_presentation == "No"))
p1 <- nrow(ed_for_log_reg %>% subset(did_the_patient_get_to_theatre_within_24_hours_of_presentation == "No" & ed == "Yes")) / n1

n2 <- nrow(ed_for_log_reg %>% subset(did_the_patient_get_to_theatre_within_24_hours_of_presentation == "Yes"))
alpha <- 0.05

p2 <- nrow(ed_for_log_reg %>% subset(did_the_patient_get_to_theatre_within_24_hours_of_presentation == "Yes" & ed == "Yes")) / n2

power <- calculate_power_table(p1, p2, n1, n2, alpha, title = "Power Calculation Results for Post-Operative Erectile Dysfunction")
power

Power Calculation Results for Post-Operative Erectile Dysfunction
Based on Two Proportion Test
Parameter	Value
Group 1 Proportion (p1)	0.22
Group 2 Proportion (p2)	0.11
Group 1 Size (n1)	76.00
Group 2 Size (n2)	138.00
Significance Level (alpha)	0.05
Calculated Power (%)	0.00

5.3.2 Logistic Regression (Firth Correction)

ed_for_log_reg <- ed_for_log_reg %>%
  left_join(knife_to_skin1, by = c("submission_id" = "submission_id")) %>%
  mutate(ed = ed %>% factor(levels = c("No", "Yes")))

ed_for_log_reg2 <- ed_for_log_reg %>%
  drop_na(
    ed,
    imaging,
    type_of_suture,
    volume_n,
    volume,
    did_the_patient_get_to_theatre_within_24_hours_of_presentation,
    what_age_was_the_patient_at_the_time_of_the_injury,
    presentation_to_knife_to_skin,
    what_age_was_the_patient_at_the_time_of_the_injury,
    was_the_patient_transferred_from_another_hospital_for_surgical_repair,
    was_a_urethral_injury_noted_at_surgery,
    laterality,
    distal,
    proximal,
    mid,
    crura,
    multi_level,
    mechanism_of_injury
  ) %>%
  mutate(
    ed = case_when(ed == "Yes" ~ 1,
                   ed == "No" ~ 0,
                   TRUE ~ NA_real_),
    presentation_to_knife_to_skin = ifelse(
    presentation_to_knife_to_skin < 0 | presentation_to_knife_to_skin > 5000,
    NA,
    presentation_to_knife_to_skin / (60*6)
  ),
    .keep = "unused"
  )

ed_lr <- logistf(
  ed ~ did_the_patient_get_to_theatre_within_24_hours_of_presentation + imaging + volume +
    what_age_was_the_patient_at_the_time_of_the_injury + was_a_urethral_injury_noted_at_surgery +
    was_the_patient_transferred_from_another_hospital_for_surgical_repair + distal + proximal + mid +
    crura + multi_level + laterality + mechanism_of_injury,
  data = ed_for_log_reg2
)

ed_lr %>%
  tbl_regression(
    exponentiate = TRUE,
    label = list(
      did_the_patient_get_to_theatre_within_24_hours_of_presentation = "Theatre within 24 hours of Presentation",
      was_a_urethral_injury_noted_at_surgery = "Urethral Injury",
      was_the_patient_transferred_from_another_hospital_for_surgical_repair = "Transfer for Surgical Repair",
      distal = "Distal Corpus Location",
    proximal = "Proximal Corpus Location",
    mid = "Mid Corpus Location",
    crura = "Crural location",
    multi_level = "Multi-level Injury",
      laterality = "Side of Injury",
      mechanism_of_injury = "Mechanism of Injury",
      what_age_was_the_patient_at_the_time_of_the_injury =
        "Age at Time of Injury",
    imaging = "Imaging",
    volume = "Centre Volume"
    )
  ) %>%
  bold_labels() %>%
  as_gt() %>%  
  tab_source_note(
    md(glue(
      "Number of patients in analysis: {nrow(ed_for_log_reg2)}"
    ))
  )

## ! `broom::tidy()` failed to tidy the model.

## ✔ `tidy_parameters()` used instead.

## ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.

Characteristic	exp(Beta)	95% CI	p-value
Theatre within 24 hours of Presentation
No	—	—
Yes	0.39	0.14, 1.03	0.056
Imaging
No	—	—
Yes	0.93	0.28, 3.23	>0.9
Centre Volume
High	—	—
Medium	1.86	0.41, 8.78	0.4
Low	1.64	0.45, 6.66	0.5
Age at Time of Injury	1.02	0.97, 1.07	0.5
Urethral Injury
No	—	—
Yes	2.14	0.63, 7.44	0.2
Transfer for Surgical Repair
No	—	—
Yes	1.53	0.46, 5.24	0.5
Distal Corpus Location
No	—	—
Yes	13.5	0.02, 37,039	0.4
Proximal Corpus Location
No	—	—
Yes	9.45	0.02, 23,178	0.5
Mid Corpus Location
No	—	—
Yes	4.33	0.01, 10,704	0.7
Crural location
No	—	—
Yes	2.60	0.01, 2,746	0.7
Multi-level Injury
No	—	—
Yes	0.64	0.00, 2,184	>0.9
Side of Injury
Bilateral	—	—
Unilateral	4.27	0.78, 47.8	0.10
Mechanism of Injury
Sexual intercourse trauma	—	—
Other	0.61	0.00, 12.2	0.8
Non-sexual blunt force trauma	0.34	0.03, 1.72	0.2
Masturbation injury	0.21	0.00, 2.01	0.2
Abbreviation: CI = Confidence Interval
Number of patients in analysis: 112

# did_the_patient_get_to_theatre_within_24_hours_of_presentation,

5.3.3 ROC

ed_lr_predict <-
  predict(ed_lr, ed_for_log_reg, type = "response")

ed_outcome <- as.factor(ifelse(ed_for_log_reg$ed == "Yes",
                                  1,
                                  0))

ed_lr_roc <-
  pROC::roc(ed_outcome, ed_lr_predict)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

auc(ed_lr_roc)

## Area under the curve: 0.7775

ci.auc(ed_lr_roc)

## 95% CI: 0.6854-0.8695 (DeLong)

ed_lr_roc_data <- cbind(sensitivities = ed_lr_roc$sensitivities,
                                specificities = ed_lr_roc$specificities) %>% as_tibble()

ggplot(as_tibble(ed_lr_roc_data), aes(x = 1 - specificities, y = sensitivities)) +
  geom_line() +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed")

5.3.4 Calibration plot

ed_lr_predict2 <- as_tibble(ed_lr_predict) %>% cbind(ed_for_log_reg) %>% subset(select = c(
  value,
  ed,
  did_the_patient_get_to_theatre_within_24_hours_of_presentation
)) %>% mutate(
  pred = value,
  .keep = "unused"
)

ed_lr_predict2$ed <- ifelse(ed_lr_predict2$ed == "Yes", 1, 0)

ed_lr_predict2$did_the_patient_get_to_theatre_within_24_hours_of_presentation <- as.integer(ed_lr_predict2$did_the_patient_get_to_theatre_within_24_hours_of_presentation)

predtools::calibration_plot(
  data = ed_lr_predict2,
  obs = "ed",
  pred = "pred",
  y_lim = c(0, 1),
  x_lim = c(0, 1),
  title = "Calibration plot for ed Model"
)

## $calibration_plot

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

5.3.5 Brier Score

brier_score <- brier(ed_lr, ed_for_log_reg2, "ed")

cbind(
  "Score" = "Brier Score",
  "Brier Score" = round(brier_score, digits = 3)
) %>% as_tibble() %>% gt() %>% gt_theme_espn()

Score	Brier Score
Brier Score	0.156

5.3.6 Nagelkerke’s R squared

pseudoR2_logistf(ed_lr, ed_for_log_reg2)

Score	Nagelkerke's R^2
Nagelkerke's R^2	0.236

5.4 Waisting / Shortening

5.4.1 Post-hoc Power calculation

n1 <- nrow(waisting_shortening_for_log_reg %>% subset(did_the_patient_get_to_theatre_within_24_hours_of_presentation == "No"))
p1 <- nrow(waisting_shortening_for_log_reg %>% subset(did_the_patient_get_to_theatre_within_24_hours_of_presentation == "No" & waisting_shortening == "Yes")) / n1

n2 <- nrow(waisting_shortening_for_log_reg %>% subset(did_the_patient_get_to_theatre_within_24_hours_of_presentation == "Yes"))
alpha <- 0.05

p2 <- nrow(waisting_shortening_for_log_reg %>% subset(did_the_patient_get_to_theatre_within_24_hours_of_presentation == "Yes" & waisting_shortening == "Yes")) / n2

power <- calculate_power_table(p1, p2, n1, n2, alpha, title = "Power Calculation Results for Post-Operative waisting_shortening")
power

Power Calculation Results for Post-Operative waisting_shortening
Based on Two Proportion Test
Parameter	Value
Group 1 Proportion (p1)	0.08
Group 2 Proportion (p2)	0.07
Group 1 Size (n1)	76.00
Group 2 Size (n2)	138.00
Significance Level (alpha)	0.05
Calculated Power (%)	1.70

5.4.2 Logistic Regression (Firth Correction)

waisting_shortening_for_log_reg <- waisting_shortening_for_log_reg %>% 
  left_join(knife_to_skin1,
            by = c("submission_id" = "submission_id")) %>%
  mutate(waisting_shortening = waisting_shortening %>% factor(levels = c("Yes",
                                           "No")))

waisting_shortening_for_log_reg2 <- waisting_shortening_for_log_reg %>%
  mutate(presentation_to_knife_to_skin = ifelse(
    presentation_to_knife_to_skin < 0 | presentation_to_knife_to_skin > 5000,
    NA,
    presentation_to_knife_to_skin / (60*6)
  ),
    waisting_shortening = case_when(waisting_shortening == "Yes" ~ 1,
                     waisting_shortening == "No" ~ 0,
                     TRUE ~ NA_real_),
  .keep = "unused"
  ) %>%
  select(
    waisting_shortening,
    imaging,
    volume_n,
    what_age_was_the_patient_at_the_time_of_the_injury,
    did_the_patient_get_to_theatre_within_24_hours_of_presentation,
    what_age_was_the_patient_at_the_time_of_the_injury,
    was_the_patient_transferred_from_another_hospital_for_surgical_repair,
    was_a_urethral_injury_noted_at_surgery,
    laterality,
    locality,
    mechanism_of_injury
  ) %>% drop_na()

waisting_shortening_lr <- logistf(
  waisting_shortening ~ did_the_patient_get_to_theatre_within_24_hours_of_presentation + imaging + volume_n +
    what_age_was_the_patient_at_the_time_of_the_injury + was_a_urethral_injury_noted_at_surgery +
    was_the_patient_transferred_from_another_hospital_for_surgical_repair + locality + laterality + mechanism_of_injury,
  data = waisting_shortening_for_log_reg2
)

waisting_shortening_lr %>%
  tbl_regression(
    exponentiate = TRUE,
    label = list(
      did_the_patient_get_to_theatre_within_24_hours_of_presentation = "Theatre within 24 hours of Presentation",
      was_a_urethral_injury_noted_at_surgery = "Urethral Injury",
      was_the_patient_transferred_from_another_hospital_for_surgical_repair = "Transfer for Surgical Repair",
      locality = "Site of Injury",
      laterality = "Side of Injury",
      mechanism_of_injury = "Mechanism of Injury",
      what_age_was_the_patient_at_the_time_of_the_injury =
        "Age at Time of Injury",
      imaging = "Imaging",
      volume_n = "Centre volume_n"
    )
  ) %>%
  bold_labels() %>%
  as_gt() %>%  
  tab_source_note(
    md(glue(
      "Number of patients in analysis: {nrow(waisting_shortening_for_log_reg2)}"
    ))
  )

## ! `broom::tidy()` failed to tidy the model.

## ✔ `tidy_parameters()` used instead.

## ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.

Characteristic	exp(Beta)	95% CI	p-value
Theatre within 24 hours of Presentation
No	—	—
Yes	0.57	0.15, 2.04	0.4
Imaging
No	—	—
Yes	0.39	0.08, 1.61	0.2
Centre volume_n	1.01	0.96, 1.06	0.7
Age at Time of Injury	1.01	0.95, 1.07	0.8
Urethral Injury
No	—	—
Yes	1.24	0.25, 4.98	0.8
Transfer for Surgical Repair
No	—	—
Yes	0.92	0.18, 4.27	>0.9
Site of Injury
Crura	—	—
Distal	0.45	0.01, 107	0.7
Mid	0.81	0.02, 173	>0.9
Multi-level	0.49	0.00, 151	0.8
Proximal	0.81	0.02, 169	>0.9
Side of Injury
Bilateral	—	—
Unilateral	0.39	0.08, 1.90	0.2
Mechanism of Injury
Sexual intercourse trauma	—	—
Other	0.85	0.01, 16.1	>0.9
Non-sexual blunt force trauma	0.27	0.00, 2.63	0.3
Masturbation injury	0.21	0.00, 3.42	0.3
Abbreviation: CI = Confidence Interval
Number of patients in analysis: 104

5.4.3 ROC

waisting_shortening_lr_predict <-
  predict(waisting_shortening_lr, waisting_shortening_for_log_reg, type = "response")

waisting_shortening_outcome <- as.factor(ifelse(waisting_shortening_for_log_reg$waisting_shortening == "Yes",
                                  1,
                                  0))

waisting_shortening_lr_roc <-
  pROC::roc(waisting_shortening_outcome, waisting_shortening_lr_predict)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

auc(waisting_shortening_lr_roc)

## Area under the curve: 0.7579

ci.auc(waisting_shortening_lr_roc)

## 95% CI: 0.6293-0.8866 (DeLong)

waisting_shortening_lr_roc_data <- cbind(sensitivities = waisting_shortening_lr_roc$sensitivities,
                                specificities = waisting_shortening_lr_roc$specificities) %>% as_tibble()

ggplot(as_tibble(waisting_shortening_lr_roc_data), aes(x = 1 - specificities, y = sensitivities)) +
  geom_line() +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed")

5.4.4 Calibration plot

waisting_shortening_lr_predict2 <- as_tibble(waisting_shortening_lr_predict) %>% cbind(waisting_shortening_for_log_reg) %>% subset(select = c(
  value,
  waisting_shortening,
  did_the_patient_get_to_theatre_within_24_hours_of_presentation
)) %>% mutate(
  pred = value,
  .keep = "unused"
)

waisting_shortening_lr_predict2$waisting_shortening <- ifelse(waisting_shortening_lr_predict2$waisting_shortening == "Yes", 1, 0)

waisting_shortening_lr_predict2$did_the_patient_get_to_theatre_within_24_hours_of_presentation <- as.integer(waisting_shortening_lr_predict2$did_the_patient_get_to_theatre_within_24_hours_of_presentation)

predtools::calibration_plot(
  data = waisting_shortening_lr_predict2,
  obs = "waisting_shortening",
  pred = "pred",
  y_lim = c(0, 1),
  x_lim = c(0, 1),
  title = "Calibration plot for waisting_shortening Model"
)

## $calibration_plot

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

5.4.5 Brier Score

brier <- function(model, data, outcome) {
  probs <- predict(model, type = "response")
  obs <- data[[outcome]]
  mean((obs - probs)^2, na.rm = TRUE)
}

brier_score <- brier(waisting_shortening_lr, waisting_shortening_for_log_reg2, "waisting_shortening") 

cbind(
  "Score" = "Brier Score",
  "Brier Score" = round(brier_score, digits = 3)
) %>% as_tibble() %>% gt() %>% gt_theme_espn()

Score	Brier Score
Brier Score	0.107

5.4.6 Nagelkerke’s R squared

pseudoR2_logistf(waisting_shortening_lr, waisting_shortening_for_log_reg2)

Score	Nagelkerke's R^2
Nagelkerke's R^2	0.146

5.5 Abnormal Curvature

5.5.1 Post-hoc Power calculation

n1 <- nrow(abnormal_curvature_for_log_reg %>% subset(did_the_patient_get_to_theatre_within_24_hours_of_presentation == "No"))
p1 <- nrow(abnormal_curvature_for_log_reg %>% subset(did_the_patient_get_to_theatre_within_24_hours_of_presentation == "No" & abnormal_curvature == "Yes")) / n1

n2 <- nrow(abnormal_curvature_for_log_reg %>% subset(did_the_patient_get_to_theatre_within_24_hours_of_presentation == "Yes"))
alpha <- 0.05

p2 <- nrow(abnormal_curvature_for_log_reg %>% subset(did_the_patient_get_to_theatre_within_24_hours_of_presentation == "Yes" & abnormal_curvature == "Yes")) / n2

power <- calculate_power_table(p1, p2, n1, n2, alpha, title = "Power Calculation Results for Post-Operative abnormal_curvature")
power

Power Calculation Results for Post-Operative abnormal_curvature
Based on Two Proportion Test
Parameter	Value
Group 1 Proportion (p1)	0.16
Group 2 Proportion (p2)	0.17
Group 1 Size (n1)	76.00
Group 2 Size (n2)	138.00
Significance Level (alpha)	0.05
Calculated Power (%)	4.90

5.5.2 Logistic Regression (Firth Correction)

Multi-level removed as model failed to separate

abnormal_curvature_for_log_reg <- abnormal_curvature_for_log_reg %>% 
  left_join(knife_to_skin1,
            by = c("submission_id" = "submission_id")) %>%
  mutate(abnormal_curvature = abnormal_curvature %>% factor(levels = c("Yes",
                                           "No")))

abnormal_curvature_for_log_reg2 <- abnormal_curvature_for_log_reg %>%
  mutate(presentation_to_knife_to_skin = ifelse(
    presentation_to_knife_to_skin < 0 | presentation_to_knife_to_skin > 5000,
    NA,
    presentation_to_knife_to_skin / (60*6)
  ),
    abnormal_curvature = case_when(abnormal_curvature == "Yes" ~ 1,
                     abnormal_curvature == "No" ~ 0,
                     TRUE ~ NA_real_),
  .keep = "unused"
  ) %>%
  select(
    abnormal_curvature,
    imaging,
    volume_n,
    what_age_was_the_patient_at_the_time_of_the_injury,
    did_the_patient_get_to_theatre_within_24_hours_of_presentation,
    what_age_was_the_patient_at_the_time_of_the_injury,
    was_the_patient_transferred_from_another_hospital_for_surgical_repair,
    was_a_urethral_injury_noted_at_surgery,
    laterality,
    distal,
    proximal,
    mid,
    crura,
    multi_level,
    mechanism_of_injury
  ) %>% drop_na()

abnormal_curvature_lr <- logistf(
  abnormal_curvature ~ did_the_patient_get_to_theatre_within_24_hours_of_presentation + imaging + volume_n + 
    what_age_was_the_patient_at_the_time_of_the_injury + was_a_urethral_injury_noted_at_surgery +
    was_the_patient_transferred_from_another_hospital_for_surgical_repair + distal + proximal + mid +
    crura + laterality + mechanism_of_injury,
  data = abnormal_curvature_for_log_reg2
)



abnormal_curvature_lr %>%
  tbl_regression(
    exponentiate = TRUE,
    label = list(
      did_the_patient_get_to_theatre_within_24_hours_of_presentation = "Theatre within 24 hours of Presentation",
      was_a_urethral_injury_noted_at_surgery = "Urethral Injury",
      was_the_patient_transferred_from_another_hospital_for_surgical_repair = "Transfer for Surgical Repair",
      locality = "Site of Injury",
      distal = "Distal Corpus Location",
    proximal = "Proximal Corpus Location",
    mid = "Mid Corpus Location",
    crura = "Crural location",
      mechanism_of_injury = "Mechanism of Injury",
      what_age_was_the_patient_at_the_time_of_the_injury =
        "Age at Time of Injury",
    imaging = "Imaging",
    volume_n = "Centre volume_n"
    )
  ) %>%
  bold_labels() %>%
  as_gt() %>%  
  tab_source_note(
    md(glue(
      "Number of patients in analysis: {nrow(abnormal_curvature_for_log_reg2)}"
    ))
  )

## ! `broom::tidy()` failed to tidy the model.

## ✔ `tidy_parameters()` used instead.

## ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.

Characteristic	exp(Beta)	95% CI	p-value
Theatre within 24 hours of Presentation
No	—	—
Yes	0.62	0.21, 1.71	0.4
Imaging
No	—	—
Yes	0.26	0.07, 0.85	0.025
Centre volume_n	1.03	0.99, 1.07	0.2
Age at Time of Injury	0.97	0.93, 1.02	0.3
Urethral Injury
No	—	—
Yes	1.79	0.57, 5.64	0.3
Transfer for Surgical Repair
No	—	—
Yes	2.49	0.79, 8.18	0.12
Distal Corpus Location
No	—	—
Yes	0.35	0.00, 5.19	0.5
Proximal Corpus Location
No	—	—
Yes	0.40	0.00, 6.88	0.6
Mid Corpus Location
No	—	—
Yes	0.22	0.00, 3.27	0.3
Crural location
No	—	—
Yes	0.10	0.00, 10.0	0.3
laterality
Bilateral	—	—
Unilateral	0.54	0.14, 2.05	0.4
Mechanism of Injury
Sexual intercourse trauma	—	—
Other	0.57	0.00, 9.04	0.7
Non-sexual blunt force trauma	0.38	0.07, 1.70	0.2
Masturbation injury	0.10	0.00, 1.16	0.069
Abbreviation: CI = Confidence Interval
Number of patients in analysis: 109

5.5.3 ROC

abnormal_curvature_lr_predict <-
  predict(abnormal_curvature_lr, abnormal_curvature_for_log_reg, type = "response")

abnormal_curvature_outcome <- as.factor(ifelse(abnormal_curvature_for_log_reg$abnormal_curvature == "Yes",
                                  1,
                                  0))

abnormal_curvature_lr_roc <-
  pROC::roc(abnormal_curvature_outcome, abnormal_curvature_lr_predict)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

auc(abnormal_curvature_lr_roc)

## Area under the curve: 0.7963

ci.auc(abnormal_curvature_lr_roc)

## 95% CI: 0.7081-0.8845 (DeLong)

abnormal_curvature_lr_roc_data <- cbind(sensitivities = abnormal_curvature_lr_roc$sensitivities,
                                specificities = abnormal_curvature_lr_roc$specificities) %>% as_tibble()

ggplot(as_tibble(abnormal_curvature_lr_roc_data), aes(x = 1 - specificities, y = sensitivities)) +
  geom_line() +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed")

5.5.4 Calibration plot

abnormal_curvature_lr_predict2 <- as_tibble(abnormal_curvature_lr_predict) %>% cbind(abnormal_curvature_for_log_reg) %>% subset(select = c(
  value,
  abnormal_curvature,
  did_the_patient_get_to_theatre_within_24_hours_of_presentation
)) %>% mutate(
  pred = value,
  .keep = "unused"
)

abnormal_curvature_lr_predict2$abnormal_curvature <- ifelse(abnormal_curvature_lr_predict2$abnormal_curvature == "Yes", 1, 0)

abnormal_curvature_lr_predict2$did_the_patient_get_to_theatre_within_24_hours_of_presentation <- as.integer(abnormal_curvature_lr_predict2$did_the_patient_get_to_theatre_within_24_hours_of_presentation)

predtools::calibration_plot(
  data = abnormal_curvature_lr_predict2,
  obs = "abnormal_curvature",
  pred = "pred",
  y_lim = c(0, 1),
  x_lim = c(0, 1),
  title = "Calibration plot for abnormal_curvature Model"
)

## $calibration_plot

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

5.5.5 Brier Score

brier <- function(model, data, outcome) {
  probs <- predict(model, type = "response")
  obs <- data[[outcome]]
  mean((obs - probs)^2, na.rm = TRUE)
}

brier_score <- brier(abnormal_curvature_lr, abnormal_curvature_for_log_reg2, "abnormal_curvature") 

cbind(
  "Score" = "Brier Score",
  "Brier Score" = round(brier_score, digits = 3)
) %>% as_tibble() %>% gt() %>% gt_theme_espn()

Score	Brier Score
Brier Score	0.165

5.5.6 Nagelkerke’s R squared

pseudoR2_logistf(abnormal_curvature_lr, abnormal_curvature_for_log_reg2)

Score	Nagelkerke's R^2
Nagelkerke's R^2	0.297

6 Logistic Regression - Time from Presentation to Theatre as Continuous Variable

6.1 Pain

6.1.1 Post-hoc Power calculation

n1 <- nrow(pain_for_log_reg %>% subset(did_the_patient_get_to_theatre_within_24_hours_of_presentation == "No"))
p1 <- nrow(pain_for_log_reg %>% subset(did_the_patient_get_to_theatre_within_24_hours_of_presentation == "No" & pain == "Yes")) / n1

n2 <- nrow(pain_for_log_reg %>% subset(did_the_patient_get_to_theatre_within_24_hours_of_presentation == "Yes"))
alpha <- 0.05

p2 <- nrow(pain_for_log_reg %>% subset(did_the_patient_get_to_theatre_within_24_hours_of_presentation == "Yes" & pain == "Yes")) / n2

power <- calculate_power_table(p1, p2, n1, n2, alpha, title = "Power Calculation Results for Post-Operative Pain")
power

Power Calculation Results for Post-Operative Pain
Based on Two Proportion Test
Parameter	Value
Group 1 Proportion (p1)	0.12
Group 2 Proportion (p2)	0.11
Group 1 Size (n1)	76.00
Group 2 Size (n2)	138.00
Significance Level (alpha)	0.05
Calculated Power (%)	1.50

6.1.2 Logistic Regression (Firth Correction)

knife_to_skin1 <- knife_to_skin %>% 
  subset(select = c(
    submission_id,
    time_from_injury_to_th
  ))

pain_for_log_reg <- pain_for_log_reg %>% 
  left_join(knife_to_skin1,
            by = c("submission_id" = "submission_id")) %>%
  mutate(pain = pain %>% factor(levels = c("Yes",
                                           "No")))

pain_for_log_reg2 <- pain_for_log_reg %>%
  mutate(presentation_to_knife_to_skin = ifelse(
    presentation_to_knife_to_skin < 0 | presentation_to_knife_to_skin > 7920,
    NA,
    presentation_to_knife_to_skin / (60*6)
  ),
    pain = case_when(pain == "Yes" ~ 1,
                     pain == "No" ~ 0,
                     TRUE ~ NA_real_),
  .keep = "unused"
  ) %>%
  select(
    pain,
    imaging,
    volume_n,
    what_age_was_the_patient_at_the_time_of_the_injury,
    presentation_to_knife_to_skin,
    what_age_was_the_patient_at_the_time_of_the_injury,
    was_the_patient_transferred_from_another_hospital_for_surgical_repair,
    was_a_urethral_injury_noted_at_surgery,
    laterality,
    distal,
    proximal,
    mid,
    crura,
    multi_level,
    mechanism_of_injury
  ) %>% drop_na()

pain_lr <- logistf(
  pain ~ presentation_to_knife_to_skin + imaging +
    what_age_was_the_patient_at_the_time_of_the_injury + was_a_urethral_injury_noted_at_surgery + volume_n +
    was_the_patient_transferred_from_another_hospital_for_surgical_repair + distal + proximal + mid +
    crura + multi_level + laterality + mechanism_of_injury,
  data = pain_for_log_reg2
)

pain_lr %>%
  tbl_regression(
    exponentiate = TRUE,
    label = list(
      presentation_to_knife_to_skin = "Time from Presentation to Theatre (6 Hour Blocks)",
      was_a_urethral_injury_noted_at_surgery = "Urethral Injury",
      was_the_patient_transferred_from_another_hospital_for_surgical_repair = "Transfer for Surgical Repair",
      laterality = "Side of Injury",
      distal = "Distal Corpus Location",
    proximal = "Proximal Corpus Location",
    mid = "Mid Corpus Location",
    crura = "Crural location",
    multi_level = "Multi-level Injury",
      mechanism_of_injury = "Mechanism of Injury",
      what_age_was_the_patient_at_the_time_of_the_injury =
        "Age at Time of Injury",
    imaging = "Imaging",
    volume_n = "Centre Volume"
    )
  ) %>%
  bold_labels() %>%
  as_gt() %>%  
  tab_source_note(
    md(glue(
      "Number of patients in analysis: {nrow(pain_for_log_reg2)}"
    ))
  )

## ! `broom::tidy()` failed to tidy the model.

## ✔ `tidy_parameters()` used instead.

## ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.

Characteristic	exp(Beta)	95% CI	p-value
Time from Presentation to Theatre (6 Hour Blocks)	0.88	0.70, 1.01	0.081
Imaging
No	—	—
Yes	1.28	0.39, 4.46	0.7
Age at Time of Injury	1.04	0.99, 1.10	0.14
Urethral Injury
No	—	—
Yes	0.86	0.18, 3.23	0.8
Centre Volume	0.93	0.88, 0.98	0.003
Transfer for Surgical Repair
No	—	—
Yes	5.32	1.57, 20.1	0.007
Distal Corpus Location
No	—	—
Yes	0.02	0.00, 94.5	0.4
Proximal Corpus Location
No	—	—
Yes	0.01	0.00, 31.1	0.3
Mid Corpus Location
No	—	—
Yes	0.01	0.00, 23.4	0.3
Crural location
No	—	—
Yes	0.21	0.00, 106	0.6
Multi-level Injury
No	—	—
Yes	408	0.04, 12,423,766	0.2
Side of Injury
Bilateral	—	—
Unilateral	0.44	0.09, 2.18	0.3
Mechanism of Injury
Sexual intercourse trauma	—	—
Other	0.74	0.00, 12.8	0.9
Non-sexual blunt force trauma	4.60	0.88, 23.4	0.069
Masturbation injury	8.83	1.07, 90.3	0.044
Abbreviation: CI = Confidence Interval
Number of patients in analysis: 125

# Add back in to get 24hr cut-off: presentation_to_knife_to_skin,
# did_the_patient_get_to_theatre_within_24_hours_of_presentation
# time_from_injury_to_th

6.1.3 ROC

pain_lr_predict <-
  predict(pain_lr, pain_for_log_reg, type = "response")

pain_outcome <- as.factor(ifelse(pain_for_log_reg$pain == "Yes",
                                  1,
                                  0))

pain_lr_roc <-
  pROC::roc(pain_outcome, pain_lr_predict)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

auc(pain_lr_roc)

## Area under the curve: 0.5595

ci.auc(pain_lr_roc)

## 95% CI: 0.4365-0.6824 (DeLong)

pain_lr_roc_data <- cbind(sensitivities = pain_lr_roc$sensitivities,
                                specificities = pain_lr_roc$specificities) %>% as_tibble()

ggplot(as_tibble(pain_lr_roc_data), aes(x = 1 - specificities, y = sensitivities)) +
  geom_line() +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed")

6.1.4 Calibration plot

pain_lr_predict2 <- as_tibble(pain_lr_predict) %>% cbind(pain_for_log_reg) %>% subset(select = c(
  value,
  pain,
  did_the_patient_get_to_theatre_within_24_hours_of_presentation
)) %>% mutate(
  pred = value,
  .keep = "unused"
)

pain_lr_predict2$pain <- ifelse(pain_lr_predict2$pain == "Yes", 1, 0)

pain_lr_predict2$did_the_patient_get_to_theatre_within_24_hours_of_presentation <- as.integer(pain_lr_predict2$did_the_patient_get_to_theatre_within_24_hours_of_presentation)

predtools::calibration_plot(
  data = pain_lr_predict2,
  obs = "pain",
  pred = "pred",
  y_lim = c(0, 1),
  x_lim = c(0, 1),
  title = "Calibration plot for Pain Model"
)

## $calibration_plot

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

6.1.5 Brier Score

brier <- function(model, data, outcome) {
  probs <- predict(model, type = "response")
  obs <- data[[outcome]]
  mean((obs - probs)^2, na.rm = TRUE)
}

brier_score <- brier(pain_lr, pain_for_log_reg2, "pain") 

cbind(
  "Score" = "Brier Score",
  "Brier Score" = round(brier_score, digits = 3)
) %>% as_tibble() %>% gt() %>% gt_theme_espn()

Score	Brier Score
Brier Score	0.106

6.1.6 Nagelkerke’s R squared

pseudoR2_logistf <- function(model, data) {
  # Extract log-likelihoods
  ll_full <- model$loglik["full"]
  ll_null <- model$loglik["null"]
  
  # Sample size
  n <- nrow(data)
  
  # Cox & Snell R²
  R2_CS <- 1 - exp((2 / n) * (ll_null - ll_full))
  
  # Nagelkerke R²
  R2_N <- R2_CS / (1 - exp((2 / n) * ll_null))

  R2_N <- R2_N %>% as_tibble()
  
  x <- cbind(
  "Score" = "Nagelkerke's R^2",
  "Nagelkerke's R^2" = round(R2_N$value, digits = 3)
) %>% as_tibble() %>% gt() %>% gt_theme_espn()
  
  # Return as a named list
  return(
    x
  )
}

pseudoR2_logistf(pain_lr, pain_for_log_reg2)

Score	Nagelkerke's R^2
Nagelkerke's R^2	0.334

6.2 ED

6.2.1 Post-hoc Power calculation

n1 <- nrow(ed_for_log_reg %>% subset(did_the_patient_get_to_theatre_within_24_hours_of_presentation == "No"))
p1 <- nrow(ed_for_log_reg %>% subset(did_the_patient_get_to_theatre_within_24_hours_of_presentation == "No" & ed == "Yes")) / n1

n2 <- nrow(ed_for_log_reg %>% subset(did_the_patient_get_to_theatre_within_24_hours_of_presentation == "Yes"))
alpha <- 0.05

p2 <- nrow(ed_for_log_reg %>% subset(did_the_patient_get_to_theatre_within_24_hours_of_presentation == "Yes" & ed == "Yes")) / n2

power <- calculate_power_table(p1, p2, n1, n2, alpha, title = "Power Calculation Results for Post-Operative Erectile Dysfunction")
power

Power Calculation Results for Post-Operative Erectile Dysfunction
Based on Two Proportion Test
Parameter	Value
Group 1 Proportion (p1)	0.22
Group 2 Proportion (p2)	0.11
Group 1 Size (n1)	76.00
Group 2 Size (n2)	138.00
Significance Level (alpha)	0.05
Calculated Power (%)	0.00

6.2.2 Logistic Regression (Firth Correction)

ed_for_log_reg <- ed_for_log_reg %>%
  left_join(knife_to_skin1, by = c("submission_id" = "submission_id")) %>%
  mutate(ed = ed %>% factor(levels = c("No", "Yes")))

ed_for_log_reg2 <- ed_for_log_reg %>%
  drop_na(
    ed,
    what_age_was_the_patient_at_the_time_of_the_injury,
    presentation_to_knife_to_skin,
    what_age_was_the_patient_at_the_time_of_the_injury,
    was_the_patient_transferred_from_another_hospital_for_surgical_repair,
    was_a_urethral_injury_noted_at_surgery,
    laterality,
    distal,
    proximal,
    mid,
    crura,
    multi_level,
    mechanism_of_injury
  ) %>%
  mutate(
    ed = case_when(ed == "Yes" ~ 1,
                   ed == "No" ~ 0,
                   TRUE ~ NA_real_),
    presentation_to_knife_to_skin = ifelse(
    presentation_to_knife_to_skin < 0 | presentation_to_knife_to_skin > 5000,
    NA,
    presentation_to_knife_to_skin / (60*6)
  ),
    .keep = "unused"
  )

ed_lr <- logistf(
  ed ~ presentation_to_knife_to_skin + imaging + volume_n +
    what_age_was_the_patient_at_the_time_of_the_injury + was_a_urethral_injury_noted_at_surgery +
    was_the_patient_transferred_from_another_hospital_for_surgical_repair + distal + proximal + mid +
    crura + multi_level + laterality + mechanism_of_injury,
  data = ed_for_log_reg2
)

ed_lr %>%
  tbl_regression(
    exponentiate = TRUE,
    label = list(
      presentation_to_knife_to_skin = "Time from Presentation to Theatre (6 Hour blocks)",
      was_a_urethral_injury_noted_at_surgery = "Urethral Injury",
      was_the_patient_transferred_from_another_hospital_for_surgical_repair = "Transfer for Surgical Repair",
      distal = "Distal Corpus Location",
    proximal = "Proximal Corpus Location",
    mid = "Mid Corpus Location",
    crura = "Crural location",
    multi_level = "Multi-level Injury",
      laterality = "Side of Injury",
      mechanism_of_injury = "Mechanism of Injury",
      what_age_was_the_patient_at_the_time_of_the_injury =
        "Age at Time of Injury",
    imaging = "Imaging",
    volume_n = "Centre Volume"
    )
  ) %>%
  bold_labels() %>%
  as_gt() %>%  
  tab_source_note(
    md(glue(
      "Number of patients in analysis: {nrow(ed_for_log_reg2)}"
    ))
  )

## ! `broom::tidy()` failed to tidy the model.

## ✔ `tidy_parameters()` used instead.

## ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.

Characteristic	exp(Beta)	95% CI	p-value
Time from Presentation to Theatre (6 Hour blocks)	1.16	0.97, 1.40	0.10
Imaging
No	—	—
Yes	0.80	0.24, 2.68	0.7
Centre Volume	0.99	0.95, 1.03	0.6
Age at Time of Injury	1.02	0.97, 1.07	0.5
Urethral Injury
No	—	—
Yes	3.18	0.94, 11.4	0.063
Transfer for Surgical Repair
No	—	—
Yes	1.41	0.42, 4.77	0.6
Distal Corpus Location
No	—	—
Yes	11.7	0.01, 38,661	0.5
Proximal Corpus Location
No	—	—
Yes	15.5	0.03, 40,330	0.4
Mid Corpus Location
No	—	—
Yes	6.13	0.01, 16,138	0.6
Crural location
No	—	—
Yes	6.42	0.02, 7,159	0.5
Multi-level Injury
No	—	—
Yes	0.20	0.00, 560	0.7
Side of Injury
Bilateral	—	—
Unilateral	7.38	1.32, 86.9	0.020
Mechanism of Injury
Sexual intercourse trauma	—	—
Other	0.42	0.00, 10.2	0.6
Non-sexual blunt force trauma	0.34	0.03, 1.71	0.2
Masturbation injury	0.37	0.00, 4.11	0.5
Abbreviation: CI = Confidence Interval
Number of patients in analysis: 136

# did_the_patient_get_to_theatre_within_24_hours_of_presentation,

6.2.3 ROC

ed_lr_predict <-
  predict(ed_lr, ed_for_log_reg, type = "response")

ed_outcome <- as.factor(ifelse(ed_for_log_reg$ed == "Yes",
                                  1,
                                  0))

ed_lr_roc <-
  pROC::roc(ed_outcome, ed_lr_predict)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

auc(ed_lr_roc)

## Area under the curve: 0.5042

ci.auc(ed_lr_roc)

## 95% CI: 0.4278-0.5805 (DeLong)

ed_lr_roc_data <- cbind(sensitivities = ed_lr_roc$sensitivities,
                                specificities = ed_lr_roc$specificities) %>% as_tibble()

ggplot(as_tibble(ed_lr_roc_data), aes(x = 1 - specificities, y = sensitivities)) +
  geom_line() +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed")

6.2.4 Calibration plot

ed_lr_predict2 <- as_tibble(ed_lr_predict) %>% cbind(ed_for_log_reg) %>% subset(select = c(
  value,
  ed,
  did_the_patient_get_to_theatre_within_24_hours_of_presentation
)) %>% mutate(
  pred = value,
  .keep = "unused"
)

ed_lr_predict2$ed <- ifelse(ed_lr_predict2$ed == "Yes", 1, 0)

ed_lr_predict2$did_the_patient_get_to_theatre_within_24_hours_of_presentation <- as.integer(ed_lr_predict2$did_the_patient_get_to_theatre_within_24_hours_of_presentation)

predtools::calibration_plot(
  data = ed_lr_predict2,
  obs = "ed",
  pred = "pred",
  y_lim = c(0, 1),
  x_lim = c(0, 1),
  title = "Calibration plot for ed Model"
)

## $calibration_plot

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

6.2.5 Brier Score

brier_score <- brier(ed_lr, ed_for_log_reg2, "ed")

## Warning in obs - probs: longer object length is not a multiple of shorter
## object length

cbind(
  "Score" = "Brier Score",
  "Brier Score" = round(brier_score, digits = 3)
) %>% as_tibble() %>% gt() %>% gt_theme_espn()

Score	Brier Score
Brier Score	0.215

6.2.6 Nagelkerke’s R squared

pseudoR2_logistf(ed_lr, ed_for_log_reg2)

Score	Nagelkerke's R^2
Nagelkerke's R^2	0.226

6.3 Waisting / Shortening

6.3.1 Post-hoc Power calculation

n1 <- nrow(waisting_shortening_for_log_reg %>% subset(did_the_patient_get_to_theatre_within_24_hours_of_presentation == "No"))
p1 <- nrow(waisting_shortening_for_log_reg %>% subset(did_the_patient_get_to_theatre_within_24_hours_of_presentation == "No" & waisting_shortening == "Yes")) / n1

n2 <- nrow(waisting_shortening_for_log_reg %>% subset(did_the_patient_get_to_theatre_within_24_hours_of_presentation == "Yes"))
alpha <- 0.05

p2 <- nrow(waisting_shortening_for_log_reg %>% subset(did_the_patient_get_to_theatre_within_24_hours_of_presentation == "Yes" & waisting_shortening == "Yes")) / n2

power <- calculate_power_table(p1, p2, n1, n2, alpha, title = "Power Calculation Results for Post-Operative waisting_shortening")
power

Power Calculation Results for Post-Operative waisting_shortening
Based on Two Proportion Test
Parameter	Value
Group 1 Proportion (p1)	0.08
Group 2 Proportion (p2)	0.07
Group 1 Size (n1)	76.00
Group 2 Size (n2)	138.00
Significance Level (alpha)	0.05
Calculated Power (%)	1.70

6.3.2 Logistic Regression (Firth Correction)

waisting_shortening_for_log_reg <- waisting_shortening_for_log_reg %>% 
  left_join(knife_to_skin1,
            by = c("submission_id" = "submission_id")) %>%
  mutate(waisting_shortening = waisting_shortening %>% factor(levels = c("Yes",
                                           "No")))

waisting_shortening_for_log_reg2 <- waisting_shortening_for_log_reg %>%
  mutate(presentation_to_knife_to_skin = ifelse(
    presentation_to_knife_to_skin < 0 | presentation_to_knife_to_skin > 5000,
    NA,
    presentation_to_knife_to_skin / (60*6)
  ),
    waisting_shortening = case_when(waisting_shortening == "Yes" ~ 1,
                     waisting_shortening == "No" ~ 0,
                     TRUE ~ NA_real_),
  .keep = "unused"
  ) %>%
  select(
    waisting_shortening,
    imaging,
    volume_n,
    what_age_was_the_patient_at_the_time_of_the_injury,
    presentation_to_knife_to_skin,
    what_age_was_the_patient_at_the_time_of_the_injury,
    was_the_patient_transferred_from_another_hospital_for_surgical_repair,
    was_a_urethral_injury_noted_at_surgery,
    laterality,
    locality,
    mechanism_of_injury
  ) %>% drop_na()

waisting_shortening_lr <- logistf(
  waisting_shortening ~ presentation_to_knife_to_skin + imaging + volume_n +
    what_age_was_the_patient_at_the_time_of_the_injury + was_a_urethral_injury_noted_at_surgery +
    was_the_patient_transferred_from_another_hospital_for_surgical_repair + locality + laterality + mechanism_of_injury,
  data = waisting_shortening_for_log_reg2
)

waisting_shortening_lr %>%
  tbl_regression(
    exponentiate = TRUE,
    label = list(
      presentation_to_knife_to_skin = "Time from Presentation to Theatre (6 Hour Blocks)",
      was_a_urethral_injury_noted_at_surgery = "Urethral Injury",
      was_the_patient_transferred_from_another_hospital_for_surgical_repair = "Transfer for Surgical Repair",
      locality = "Site of Injury",
      laterality = "Side of Injury",
      mechanism_of_injury = "Mechanism of Injury",
      what_age_was_the_patient_at_the_time_of_the_injury =
        "Age at Time of Injury",
      imaging = "Imaging",
      volume_n = "Centre Volume"
    )
  ) %>%
  bold_labels() %>%
  as_gt() %>%  
  tab_source_note(
    md(glue(
      "Number of patients in analysis: {nrow(waisting_shortening_for_log_reg2)}"
    ))
  )

## ! `broom::tidy()` failed to tidy the model.

## ✔ `tidy_parameters()` used instead.

## ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.

Characteristic	exp(Beta)	95% CI	p-value
Time from Presentation to Theatre (6 Hour Blocks)	0.95	0.69, 1.26	0.7
Imaging
No	—	—
Yes	0.47	0.08, 2.26	0.3
Centre Volume	1.02	0.97, 1.08	0.4
Age at Time of Injury	1.02	0.96, 1.09	0.6
Urethral Injury
No	—	—
Yes	1.72	0.31, 7.35	0.5
Transfer for Surgical Repair
No	—	—
Yes	0.59	0.09, 3.04	0.5
Site of Injury
Crura	—	—
Distal	0.10	0.00, 30.2	0.4
Mid	0.36	0.01, 76.1	0.6
Multi-level	0.36	0.00, 106	0.7
Proximal	0.57	0.01, 116	0.8
Side of Injury
Bilateral	—	—
Unilateral	0.25	0.05, 1.14	0.072
Mechanism of Injury
Sexual intercourse trauma	—	—
Other	0.45	0.00, 8.59	0.6
Non-sexual blunt force trauma	0.30	0.00, 3.13	0.4
Masturbation injury	2.03	0.01, 41.8	0.7
Abbreviation: CI = Confidence Interval
Number of patients in analysis: 101

6.3.3 ROC

waisting_shortening_lr_predict <-
  predict(waisting_shortening_lr, waisting_shortening_for_log_reg, type = "response")

waisting_shortening_outcome <- as.factor(ifelse(waisting_shortening_for_log_reg$waisting_shortening == "Yes",
                                  1,
                                  0))

waisting_shortening_lr_roc <-
  pROC::roc(waisting_shortening_outcome, waisting_shortening_lr_predict)

## Setting levels: control = 0, case = 1

## Setting direction: controls > cases

auc(waisting_shortening_lr_roc)

## Area under the curve: 0.5429

ci.auc(waisting_shortening_lr_roc)

## 95% CI: 0.3802-0.7057 (DeLong)

waisting_shortening_lr_roc_data <- cbind(sensitivities = waisting_shortening_lr_roc$sensitivities,
                                specificities = waisting_shortening_lr_roc$specificities) %>% as_tibble()

ggplot(as_tibble(waisting_shortening_lr_roc_data), aes(x = 1 - specificities, y = sensitivities)) +
  geom_line() +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed")

6.3.4 Calibration plot

waisting_shortening_lr_predict2 <- as_tibble(waisting_shortening_lr_predict) %>% cbind(waisting_shortening_for_log_reg) %>% subset(select = c(
  value,
  waisting_shortening,
  did_the_patient_get_to_theatre_within_24_hours_of_presentation
)) %>% mutate(
  pred = value,
  .keep = "unused"
)

waisting_shortening_lr_predict2$waisting_shortening <- ifelse(waisting_shortening_lr_predict2$waisting_shortening == "Yes", 1, 0)

waisting_shortening_lr_predict2$did_the_patient_get_to_theatre_within_24_hours_of_presentation <- as.integer(waisting_shortening_lr_predict2$did_the_patient_get_to_theatre_within_24_hours_of_presentation)

predtools::calibration_plot(
  data = waisting_shortening_lr_predict2,
  obs = "waisting_shortening",
  pred = "pred",
  y_lim = c(0, 1),
  x_lim = c(0, 1),
  title = "Calibration plot for waisting_shortening Model"
)

## $calibration_plot

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

6.3.5 Brier Score

brier <- function(model, data, outcome) {
  probs <- predict(model, type = "response")
  obs <- data[[outcome]]
  mean((obs - probs)^2, na.rm = TRUE)
}

brier_score <- brier(waisting_shortening_lr, waisting_shortening_for_log_reg2, "waisting_shortening") 

cbind(
  "Score" = "Brier Score",
  "Brier Score" = round(brier_score, digits = 3)
) %>% as_tibble() %>% gt() %>% gt_theme_espn()

Score	Brier Score
Brier Score	0.094

6.3.6 Nagelkerke’s R squared

pseudoR2_logistf(waisting_shortening_lr, waisting_shortening_for_log_reg2)

Score	Nagelkerke's R^2
Nagelkerke's R^2	0.2

6.4 Abnormal Curvature

6.4.1 Post-hoc Power calculation

n1 <- nrow(abnormal_curvature_for_log_reg %>% subset(did_the_patient_get_to_theatre_within_24_hours_of_presentation == "No"))
p1 <- nrow(abnormal_curvature_for_log_reg %>% subset(did_the_patient_get_to_theatre_within_24_hours_of_presentation == "No" & abnormal_curvature == "Yes")) / n1

n2 <- nrow(abnormal_curvature_for_log_reg %>% subset(did_the_patient_get_to_theatre_within_24_hours_of_presentation == "Yes"))
alpha <- 0.05

p2 <- nrow(abnormal_curvature_for_log_reg %>% subset(did_the_patient_get_to_theatre_within_24_hours_of_presentation == "Yes" & abnormal_curvature == "Yes")) / n2

power <- calculate_power_table(p1, p2, n1, n2, alpha, title = "Power Calculation Results for Post-Operative abnormal_curvature")
power

Power Calculation Results for Post-Operative abnormal_curvature
Based on Two Proportion Test
Parameter	Value
Group 1 Proportion (p1)	0.16
Group 2 Proportion (p2)	0.17
Group 1 Size (n1)	76.00
Group 2 Size (n2)	138.00
Significance Level (alpha)	0.05
Calculated Power (%)	4.90

6.4.2 Logistic Regression (Firth Correction)

Multi-level removed as model would not separate

abnormal_curvature_for_log_reg <- abnormal_curvature_for_log_reg %>% 
  left_join(knife_to_skin1,
            by = c("submission_id" = "submission_id")) %>%
  mutate(abnormal_curvature = abnormal_curvature %>% factor(levels = c("Yes",
                                           "No")))

abnormal_curvature_for_log_reg2 <- abnormal_curvature_for_log_reg %>%
  mutate(presentation_to_knife_to_skin = ifelse(
    presentation_to_knife_to_skin < 0 | presentation_to_knife_to_skin > 5000,
    NA,
    presentation_to_knife_to_skin / (60*6)
  ),
    abnormal_curvature = case_when(abnormal_curvature == "Yes" ~ 1,
                     abnormal_curvature == "No" ~ 0,
                     TRUE ~ NA_real_),
  .keep = "unused"
  ) %>%
  select(
    abnormal_curvature,
    imaging,
    volume_n,
    what_age_was_the_patient_at_the_time_of_the_injury,
    presentation_to_knife_to_skin,
    what_age_was_the_patient_at_the_time_of_the_injury,
    was_the_patient_transferred_from_another_hospital_for_surgical_repair,
    was_a_urethral_injury_noted_at_surgery,
    laterality,
    distal,
    proximal,
    mid,
    crura,
    multi_level,
    mechanism_of_injury
  ) %>% drop_na()

abnormal_curvature_lr <- logistf(
  abnormal_curvature ~ presentation_to_knife_to_skin + imaging + volume_n +
    what_age_was_the_patient_at_the_time_of_the_injury + was_a_urethral_injury_noted_at_surgery +
    was_the_patient_transferred_from_another_hospital_for_surgical_repair + distal + proximal + mid +
    crura + laterality + mechanism_of_injury,
  data = abnormal_curvature_for_log_reg2
)

abnormal_curvature_lr %>%
  tbl_regression(
    exponentiate = TRUE,
    label = list(
      presentation_to_knife_to_skin = "Time from Presentation to Theatre (6 Hour Blocks)",
      was_a_urethral_injury_noted_at_surgery = "Urethral Injury",
      was_the_patient_transferred_from_another_hospital_for_surgical_repair = "Transfer for Surgical Repair",
      laterality = "Side of Injury",
      distal = "Distal Corpus Location",
    proximal = "Proximal Corpus Location",
    mid = "Mid Corpus Location",
    crura = "Crural location",
    multi_level = "Multi-level Injury",
      mechanism_of_injury = "Mechanism of Injury",
      what_age_was_the_patient_at_the_time_of_the_injury =
        "Age at Time of Injury",
    imaging = "Imaging",
    volume_n = "Centre Volume"
    )
  ) %>%
  bold_labels() %>%
  as_gt() %>%  
  tab_source_note(
    md(glue(
      "Number of patients in analysis: {nrow(abnormal_curvature_for_log_reg2)}"
    ))
  )

## ! `broom::tidy()` failed to tidy the model.

## ✔ `tidy_parameters()` used instead.

## ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.

Characteristic	exp(Beta)	95% CI	p-value
Time from Presentation to Theatre (6 Hour Blocks)	1.05	0.88, 1.26	0.6
Imaging
No	—	—
Yes	0.37	0.11, 1.19	0.10
Centre Volume	1.02	0.98, 1.07	0.2
Age at Time of Injury	0.97	0.93, 1.02	0.2
Urethral Injury
No	—	—
Yes	1.11	0.33, 3.52	0.9
Transfer for Surgical Repair
No	—	—
Yes	2.14	0.71, 6.66	0.2
Distal Corpus Location
No	—	—
Yes	0.41	0.00, 5.88	0.6
Proximal Corpus Location
No	—	—
Yes	0.32	0.00, 5.43	0.5
Mid Corpus Location
No	—	—
Yes	0.19	0.00, 2.76	0.2
Crural location
No	—	—
Yes	0.09	0.00, 8.13	0.3
Side of Injury
Bilateral	—	—
Unilateral	0.52	0.14, 1.84	0.3
Mechanism of Injury
Sexual intercourse trauma	—	—
Other	0.56	0.00, 8.06	0.7
Non-sexual blunt force trauma	0.35	0.06, 1.54	0.2
Masturbation injury	0.20	0.00, 2.32	0.2
Abbreviation: CI = Confidence Interval
Number of patients in analysis: 107

6.4.3 ROC

abnormal_curvature_lr_predict <-
  predict(abnormal_curvature_lr, abnormal_curvature_for_log_reg, type = "response")

abnormal_curvature_outcome <- as.factor(ifelse(abnormal_curvature_for_log_reg$abnormal_curvature == "Yes",
                                  1,
                                  0))

abnormal_curvature_lr_roc <-
  pROC::roc(abnormal_curvature_outcome, abnormal_curvature_lr_predict)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

auc(abnormal_curvature_lr_roc)

## Area under the curve: 0.5337

ci.auc(abnormal_curvature_lr_roc)

## 95% CI: 0.4327-0.6347 (DeLong)

abnormal_curvature_lr_roc_data <- cbind(sensitivities = abnormal_curvature_lr_roc$sensitivities,
                                specificities = abnormal_curvature_lr_roc$specificities) %>% as_tibble()

ggplot(as_tibble(abnormal_curvature_lr_roc_data), aes(x = 1 - specificities, y = sensitivities)) +
  geom_line() +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed")

6.4.4 Calibration plot

abnormal_curvature_lr_predict2 <- as_tibble(abnormal_curvature_lr_predict) %>% cbind(abnormal_curvature_for_log_reg) %>% subset(select = c(
  value,
  abnormal_curvature,
  did_the_patient_get_to_theatre_within_24_hours_of_presentation
)) %>% mutate(
  pred = value,
  .keep = "unused"
)

abnormal_curvature_lr_predict2$abnormal_curvature <- ifelse(abnormal_curvature_lr_predict2$abnormal_curvature == "Yes", 1, 0)

abnormal_curvature_lr_predict2$did_the_patient_get_to_theatre_within_24_hours_of_presentation <- as.integer(abnormal_curvature_lr_predict2$did_the_patient_get_to_theatre_within_24_hours_of_presentation)

predtools::calibration_plot(
  data = abnormal_curvature_lr_predict2,
  obs = "abnormal_curvature",
  pred = "pred",
  y_lim = c(0, 1),
  x_lim = c(0, 1),
  title = "Calibration plot for abnormal_curvature Model"
)

## $calibration_plot

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

6.4.5 Brier Score

brier <- function(model, data, outcome) {
  probs <- predict(model, type = "response")
  obs <- data[[outcome]]
  mean((obs - probs)^2, na.rm = TRUE)
}

brier_score <- brier(abnormal_curvature_lr, abnormal_curvature_for_log_reg2, "abnormal_curvature") 

cbind(
  "Score" = "Brier Score",
  "Brier Score" = round(brier_score, digits = 3)
) %>% as_tibble() %>% gt() %>% gt_theme_espn()

Score	Brier Score
Brier Score	0.174

6.4.6 Nagelkerke’s R squared

pseudoR2_logistf(abnormal_curvature_lr, abnormal_curvature_for_log_reg2)

Score	Nagelkerke's R^2
Nagelkerke's R^2	0.266

7 Logistic Regression - Time from Presentation to Theatre with 18 hour cut-off

7.1 Pain

7.1.1 Post-hoc Power calculation

n1 <- nrow(pain_for_log_reg_18 %>% subset(did_the_patient_get_to_theatre_within_18_hours_of_presentation == "No"))
p1 <- nrow(pain_for_log_reg_18 %>% subset(did_the_patient_get_to_theatre_within_18_hours_of_presentation == "No" & pain == "Yes")) / n1

n2 <- nrow(pain_for_log_reg_18 %>% subset(did_the_patient_get_to_theatre_within_18_hours_of_presentation == "Yes"))
alpha <- 0.05

p2 <- nrow(pain_for_log_reg_18 %>% subset(did_the_patient_get_to_theatre_within_18_hours_of_presentation == "Yes" & pain == "Yes")) / n2

power <- calculate_power_table(p1, p2, n1, n2, alpha, title = "Power Calculation Results for Post-Operative Pain")
power

Power Calculation Results for Post-Operative Pain
Based on Two Proportion Test
Parameter	Value
Group 1 Proportion (p1)	0.11
Group 2 Proportion (p2)	0.12
Group 1 Size (n1)	112.00
Group 2 Size (n2)	133.00
Significance Level (alpha)	0.05
Calculated Power (%)	5.10

7.1.2 Logistic Regression (Firth Correction)

knife_to_skin1 <- knife_to_skin %>% 
  subset(select = c(
    submission_id,
    time_from_injury_to_th
  ))

pain_for_log_reg_18 <- pain_for_log_reg_18 %>% 
  left_join(knife_to_skin1,
            by = c("submission_id" = "submission_id")) %>%
  mutate(pain = pain %>% factor(levels = c("Yes",
                                           "No")))

pain_for_log_reg_182 <- pain_for_log_reg_18 %>%
  mutate(presentation_to_knife_to_skin = ifelse(
    presentation_to_knife_to_skin < 0 | presentation_to_knife_to_skin > 7920,
    NA,
    presentation_to_knife_to_skin / (60*6)
  ),
    pain = case_when(pain == "Yes" ~ 1,
                     pain == "No" ~ 0,
                     TRUE ~ NA_real_),
  .keep = "unused"
  ) %>%
  select(
    pain,
    imaging,
    volume_n,
    what_age_was_the_patient_at_the_time_of_the_injury,
    did_the_patient_get_to_theatre_within_18_hours_of_presentation,
    what_age_was_the_patient_at_the_time_of_the_injury,
    was_the_patient_transferred_from_another_hospital_for_surgical_repair,
    was_a_urethral_injury_noted_at_surgery,
    laterality,
    distal,
    proximal,
    mid,
    crura,
    multi_level,
    mechanism_of_injury
  ) %>% drop_na()

pain_lr <- logistf(
  pain ~ did_the_patient_get_to_theatre_within_18_hours_of_presentation + imaging + volume_n +
    what_age_was_the_patient_at_the_time_of_the_injury + was_a_urethral_injury_noted_at_surgery +
    was_the_patient_transferred_from_another_hospital_for_surgical_repair + distal + proximal + mid +
    crura + multi_level + laterality + mechanism_of_injury,
  data = pain_for_log_reg_182
)

pain_lr %>%
  tbl_regression(
    exponentiate = TRUE,
    label = list(
      did_the_patient_get_to_theatre_within_18_hours_of_presentation = "Theatre within 18 hours of Presentation",
      was_a_urethral_injury_noted_at_surgery = "Urethral Injury",
      was_the_patient_transferred_from_another_hospital_for_surgical_repair = "Transfer for Surgical Repair",
      laterality = "Side of Injury",
      distal = "Distal Corpus Location",
    proximal = "Proximal Corpus Location",
    mid = "Mid Corpus Location",
    crura = "Crural location",
    multi_level = "Multi-level Injury",
      mechanism_of_injury = "Mechanism of Injury",
      what_age_was_the_patient_at_the_time_of_the_injury =
        "Age at Time of Injury",
    imaging = "Imaging",
    volume_n = "Centre Volume"
    )
  ) %>%
  bold_labels() %>%
  as_gt() %>%  
  tab_source_note(
    md(glue(
      "Number of patients in analysis: {nrow(pain_for_log_reg_182)}"
    ))
  )

## ! `broom::tidy()` failed to tidy the model.

## ✔ `tidy_parameters()` used instead.

## ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.

Characteristic	exp(Beta)	95% CI	p-value
Theatre within 18 hours of Presentation
No	—	—
Yes	1.51	0.51, 4.69	0.5
Imaging
No	—	—
Yes	1.05	0.36, 3.18	>0.9
Centre Volume	0.94	0.89, 0.99	0.009
Age at Time of Injury	1.03	0.99, 1.08	0.2
Urethral Injury
No	—	—
Yes	1.16	0.32, 3.73	0.8
Transfer for Surgical Repair
No	—	—
Yes	4.26	1.32, 14.8	0.015
Distal Corpus Location
No	—	—
Yes	0.00	0.00, 7.11	0.14
Proximal Corpus Location
No	—	—
Yes	0.01	0.00, 8.13	0.2
Mid Corpus Location
No	—	—
Yes	0.00	0.00, 4.75	0.12
Crural location
No	—	—
Yes	0.07	0.00, 32.8	0.4
Multi-level Injury
No	—	—
Yes	1,429	0.25, 20,142,183	0.10
Side of Injury
Bilateral	—	—
Unilateral	0.75	0.19, 3.45	0.7
Mechanism of Injury
Sexual intercourse trauma	—	—
Other	0.97	0.01, 12.1	>0.9
Non-sexual blunt force trauma	3.91	0.79, 17.5	0.091
Masturbation injury	4.89	0.65, 37.6	0.12
Abbreviation: CI = Confidence Interval
Number of patients in analysis: 141

# Add back in to get 24hr cut-off: presentation_to_knife_to_skin,
# did_the_patient_get_to_theatre_within_18_hours_of_presentation
# time_from_injury_to_th

7.1.3 ROC

pain_lr_predict <-
  predict(pain_lr, pain_for_log_reg_18, type = "response")

pain_outcome <- as.factor(ifelse(pain_for_log_reg_18$pain == "Yes",
                                  1,
                                  0))

pain_lr_roc <-
  pROC::roc(pain_outcome, pain_lr_predict)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

auc(pain_lr_roc)

## Area under the curve: 0.7865

ci.auc(pain_lr_roc)

## 95% CI: 0.6891-0.8839 (DeLong)

pain_lr_roc_data <- cbind(sensitivities = pain_lr_roc$sensitivities,
                                specificities = pain_lr_roc$specificities) %>% as_tibble()

ggplot(as_tibble(pain_lr_roc_data), aes(x = 1 - specificities, y = sensitivities)) +
  geom_line() +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed")

7.1.4 Calibration plot

pain_lr_predict2 <- as_tibble(pain_lr_predict) %>% cbind(pain_for_log_reg_18) %>% subset(select = c(
  value,
  pain,
  did_the_patient_get_to_theatre_within_18_hours_of_presentation
)) %>% mutate(
  pred = value,
  .keep = "unused"
)

pain_lr_predict2$pain <- ifelse(pain_lr_predict2$pain == "Yes", 1, 0)

pain_lr_predict2$did_the_patient_get_to_theatre_within_18_hours_of_presentation <- as.integer(pain_lr_predict2$did_the_patient_get_to_theatre_within_18_hours_of_presentation)

## Warning: NAs introduced by coercion

predtools::calibration_plot(
  data = pain_lr_predict2,
  obs = "pain",
  pred = "pred",
  y_lim = c(0, 1),
  x_lim = c(0, 1),
  title = "Calibration plot for Pain Model"
)

## $calibration_plot

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

7.1.5 Brier Score

brier <- function(model, data, outcome) {
  probs <- predict(model, type = "response")
  obs <- data[[outcome]]
  mean((obs - probs)^2, na.rm = TRUE)
}

brier_score <- brier(pain_lr, pain_for_log_reg_182, "pain") 

cbind(
  "Score" = "Brier Score",
  "Brier Score" = round(brier_score, digits = 3)
) %>% as_tibble() %>% gt() %>% gt_theme_espn()

Score	Brier Score
Brier Score	0.118

7.1.6 Nagelkerke’s R squared

pseudoR2_logistf <- function(model, data) {
  # Extract log-likelihoods
  ll_full <- model$loglik["full"]
  ll_null <- model$loglik["null"]
  
  # Sample size
  n <- nrow(data)
  
  # Cox & Snell R²
  R2_CS <- 1 - exp((2 / n) * (ll_null - ll_full))
  
  # Nagelkerke R²
  R2_N <- R2_CS / (1 - exp((2 / n) * ll_null))

  R2_N <- R2_N %>% as_tibble()
  
  x <- cbind(
  "Score" = "Nagelkerke's R^2",
  "Nagelkerke's R^2" = round(R2_N$value, digits = 3)
) %>% as_tibble() %>% gt() %>% gt_theme_espn()
  
  # Return as a named list
  return(
    x
  )
}

pseudoR2_logistf(pain_lr, pain_for_log_reg_182)

Score	Nagelkerke's R^2
Nagelkerke's R^2	0.25

7.2 ED

7.2.1 Post-hoc Power calculation

n1 <- nrow(ed_for_log_reg_18 %>% subset(did_the_patient_get_to_theatre_within_18_hours_of_presentation == "No"))
p1 <- nrow(ed_for_log_reg_18 %>% subset(did_the_patient_get_to_theatre_within_18_hours_of_presentation == "No" & ed == "Yes")) / n1

n2 <- nrow(ed_for_log_reg_18 %>% subset(did_the_patient_get_to_theatre_within_18_hours_of_presentation == "Yes"))
alpha <- 0.05

p2 <- nrow(ed_for_log_reg_18 %>% subset(did_the_patient_get_to_theatre_within_18_hours_of_presentation == "Yes" & ed == "Yes")) / n2

power <- calculate_power_table(p1, p2, n1, n2, alpha, title = "Power Calculation Results for Post-Operative Erectile Dysfunction")
power

Power Calculation Results for Post-Operative Erectile Dysfunction
Based on Two Proportion Test
Parameter	Value
Group 1 Proportion (p1)	0.16
Group 2 Proportion (p2)	0.16
Group 1 Size (n1)	112.00
Group 2 Size (n2)	133.00
Significance Level (alpha)	0.05
Calculated Power (%)	2.20

7.2.2 Logistic Regression (Firth Correction)

ed_for_log_reg_18 <- ed_for_log_reg_18 %>%
  left_join(knife_to_skin1, by = c("submission_id" = "submission_id")) %>%
  mutate(ed = ed %>% factor(levels = c("No", "Yes")))

ed_for_log_reg_182 <- ed_for_log_reg_18 %>%
  drop_na(
    ed,
    imaging,
    volume_n,
    did_the_patient_get_to_theatre_within_18_hours_of_presentation,
    what_age_was_the_patient_at_the_time_of_the_injury,
    presentation_to_knife_to_skin,
    what_age_was_the_patient_at_the_time_of_the_injury,
    was_the_patient_transferred_from_another_hospital_for_surgical_repair,
    was_a_urethral_injury_noted_at_surgery,
    laterality,
    distal,
    proximal,
    mid,
    crura,
    multi_level,
    mechanism_of_injury
  ) %>%
  mutate(
    ed = case_when(ed == "Yes" ~ 1,
                   ed == "No" ~ 0,
                   TRUE ~ NA_real_),
    presentation_to_knife_to_skin = ifelse(
    presentation_to_knife_to_skin < 0 | presentation_to_knife_to_skin > 5000,
    NA,
    presentation_to_knife_to_skin / (60*6)
  ),
    .keep = "unused"
  )

ed_lr <- logistf(
  ed ~ did_the_patient_get_to_theatre_within_18_hours_of_presentation + imaging + volume_n +
    what_age_was_the_patient_at_the_time_of_the_injury + was_a_urethral_injury_noted_at_surgery +
    was_the_patient_transferred_from_another_hospital_for_surgical_repair + distal + proximal + mid +
    crura + multi_level + laterality + mechanism_of_injury,
  data = ed_for_log_reg_182
)

ed_lr %>%
  tbl_regression(
    exponentiate = TRUE,
    label = list(
      did_the_patient_get_to_theatre_within_18_hours_of_presentation = "Theatre within 18 hours of Presentation",
      was_a_urethral_injury_noted_at_surgery = "Urethral Injury",
      was_the_patient_transferred_from_another_hospital_for_surgical_repair = "Transfer for Surgical Repair",
      distal = "Distal Corpus Location",
    proximal = "Proximal Corpus Location",
    mid = "Mid Corpus Location",
    crura = "Crural location",
    multi_level = "Multi-level Injury",
      laterality = "Side of Injury",
      mechanism_of_injury = "Mechanism of Injury",
      what_age_was_the_patient_at_the_time_of_the_injury =
        "Age at Time of Injury",
    imaging = "Imaging",
    volume_n = "Centre Volume"
    )
  ) %>%
  bold_labels() %>%
  as_gt() %>%  
  tab_source_note(
    md(glue(
      "Number of patients in analysis: {nrow(ed_for_log_reg_182)}"
    ))
  )

## ! `broom::tidy()` failed to tidy the model.

## ✔ `tidy_parameters()` used instead.

## ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.

Characteristic	exp(Beta)	95% CI	p-value
Theatre within 18 hours of Presentation
No	—	—
Yes	0.78	0.31, 1.92	0.6
Imaging
No	—	—
Yes	1.36	0.48, 3.97	0.6
Centre Volume	0.98	0.94, 1.01	0.3
Age at Time of Injury	1.02	0.98, 1.07	0.3
Urethral Injury
No	—	—
Yes	2.64	0.91, 7.93	0.074
Transfer for Surgical Repair
No	—	—
Yes	1.76	0.56, 5.83	0.3
Distal Corpus Location
No	—	—
Yes	99.6	0.17, 256,830	0.2
Proximal Corpus Location
No	—	—
Yes	63.6	0.15, 130,045	0.2
Mid Corpus Location
No	—	—
Yes	30.0	0.07, 61,226	0.3
Crural location
No	—	—
Yes	21.3	0.11, 19,175	0.3
Multi-level Injury
No	—	—
Yes	0.03	0.00, 55.2	0.4
Side of Injury
Bilateral	—	—
Unilateral	5.10	1.16, 33.3	0.029
Mechanism of Injury
Sexual intercourse trauma	—	—
Other	0.28	0.00, 3.88	0.4
Non-sexual blunt force trauma	0.25	0.03, 1.20	0.087
Masturbation injury	0.14	0.00, 1.27	0.088
Abbreviation: CI = Confidence Interval
Number of patients in analysis: 136

# did_the_patient_get_to_theatre_within_18_hours_of_presentation,

7.2.3 ROC

ed_lr_predict <-
  predict(ed_lr, ed_for_log_reg_18, type = "response")

ed_outcome <- as.factor(ifelse(ed_for_log_reg_18$ed == "Yes",
                                  1,
                                  0))

ed_lr_roc <-
  pROC::roc(ed_outcome, ed_lr_predict)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

auc(ed_lr_roc)

## Area under the curve: 0.7903

ci.auc(ed_lr_roc)

## 95% CI: 0.7089-0.8717 (DeLong)

ed_lr_roc_data <- cbind(sensitivities = ed_lr_roc$sensitivities,
                                specificities = ed_lr_roc$specificities) %>% as_tibble()

ggplot(as_tibble(ed_lr_roc_data), aes(x = 1 - specificities, y = sensitivities)) +
  geom_line() +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed")

7.2.4 Calibration plot

ed_lr_predict2 <- as_tibble(ed_lr_predict) %>% cbind(ed_for_log_reg_18) %>% subset(select = c(
  value,
  ed,
  did_the_patient_get_to_theatre_within_18_hours_of_presentation
)) %>% mutate(
  pred = value,
  .keep = "unused"
)

ed_lr_predict2$ed <- ifelse(ed_lr_predict2$ed == "Yes", 1, 0)

ed_lr_predict2$did_the_patient_get_to_theatre_within_18_hours_of_presentation <- as.integer(ed_lr_predict2$did_the_patient_get_to_theatre_within_18_hours_of_presentation)

## Warning: NAs introduced by coercion

predtools::calibration_plot(
  data = ed_lr_predict2,
  obs = "ed",
  pred = "pred",
  y_lim = c(0, 1),
  x_lim = c(0, 1),
  title = "Calibration plot for ed Model"
)

## $calibration_plot

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

7.2.5 Brier Score

brier_score <- brier(ed_lr, ed_for_log_reg_182, "ed")

cbind(
  "Score" = "Brier Score",
  "Brier Score" = round(brier_score, digits = 3)
) %>% as_tibble() %>% gt() %>% gt_theme_espn()

Score	Brier Score
Brier Score	0.161

7.2.6 Nagelkerke’s R squared

pseudoR2_logistf(ed_lr, ed_for_log_reg_182)

Score	Nagelkerke's R^2
Nagelkerke's R^2	0.232

7.3 Waisting / Shortening

7.3.1 Post-hoc Power calculation

n1 <- nrow(waisting_shortening_for_log_reg_18 %>% subset(did_the_patient_get_to_theatre_within_18_hours_of_presentation == "No"))
p1 <- nrow(waisting_shortening_for_log_reg_18 %>% subset(did_the_patient_get_to_theatre_within_18_hours_of_presentation == "No" & waisting_shortening == "Yes")) / n1

n2 <- nrow(waisting_shortening_for_log_reg_18 %>% subset(did_the_patient_get_to_theatre_within_18_hours_of_presentation == "Yes"))
alpha <- 0.05

p2 <- nrow(waisting_shortening_for_log_reg_18 %>% subset(did_the_patient_get_to_theatre_within_18_hours_of_presentation == "Yes" & waisting_shortening == "Yes")) / n2

power <- calculate_power_table(p1, p2, n1, n2, alpha, title = "Power Calculation Results for Post-Operative waisting_shortening")
power

Power Calculation Results for Post-Operative waisting_shortening
Based on Two Proportion Test
Parameter	Value
Group 1 Proportion (p1)	0.09
Group 2 Proportion (p2)	0.05
Group 1 Size (n1)	112.00
Group 2 Size (n2)	133.00
Significance Level (alpha)	0.05
Calculated Power (%)	0.10

7.3.2 Logistic Regression (Firth Correction)

waisting_shortening_for_log_reg_18 <- waisting_shortening_for_log_reg_18 %>% 
  left_join(knife_to_skin1,
            by = c("submission_id" = "submission_id")) %>%
  mutate(waisting_shortening = waisting_shortening %>% factor(levels = c("Yes",
                                           "No")))

waisting_shortening_for_log_reg_182 <- waisting_shortening_for_log_reg_18 %>%
  mutate(presentation_to_knife_to_skin = ifelse(
    presentation_to_knife_to_skin < 0 | presentation_to_knife_to_skin > 5000,
    NA,
    presentation_to_knife_to_skin / (60*6)
  ),
    waisting_shortening = case_when(waisting_shortening == "Yes" ~ 1,
                     waisting_shortening == "No" ~ 0,
                     TRUE ~ NA_real_),
  .keep = "unused"
  ) %>%
  select(
    waisting_shortening,
    imaging,
    volume_n,
    what_age_was_the_patient_at_the_time_of_the_injury,
    did_the_patient_get_to_theatre_within_18_hours_of_presentation,
    what_age_was_the_patient_at_the_time_of_the_injury,
    was_the_patient_transferred_from_another_hospital_for_surgical_repair,
    was_a_urethral_injury_noted_at_surgery,
    laterality,
    locality,
    mechanism_of_injury
  ) %>% drop_na()

waisting_shortening_lr <- logistf(
  waisting_shortening ~ did_the_patient_get_to_theatre_within_18_hours_of_presentation + imaging + volume_n +
    what_age_was_the_patient_at_the_time_of_the_injury + was_a_urethral_injury_noted_at_surgery +
    was_the_patient_transferred_from_another_hospital_for_surgical_repair + locality + laterality + mechanism_of_injury,
  data = waisting_shortening_for_log_reg_182
)

waisting_shortening_lr %>%
  tbl_regression(
    exponentiate = TRUE,
    label = list(
      did_the_patient_get_to_theatre_within_18_hours_of_presentation = "Theatre within 18 hours of Presentation",
      was_a_urethral_injury_noted_at_surgery = "Urethral Injury",
      was_the_patient_transferred_from_another_hospital_for_surgical_repair = "Transfer for Surgical Repair",
      locality = "Site of Injury",
      laterality = "Side of Injury",
      mechanism_of_injury = "Mechanism of Injury",
      what_age_was_the_patient_at_the_time_of_the_injury =
        "Age at Time of Injury",
      imaging = "Imaging",
      volume_n = "Centre Volume"
    )
  ) %>%
  bold_labels() %>%
  as_gt() %>%  
  tab_source_note(
    md(glue(
      "Number of patients in analysis: {nrow(waisting_shortening_for_log_reg_182)}"
    ))
  )

## ! `broom::tidy()` failed to tidy the model.

## ✔ `tidy_parameters()` used instead.

## ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.

Characteristic	exp(Beta)	95% CI	p-value
Theatre within 18 hours of Presentation
No	—	—
Yes	0.44	0.11, 1.51	0.2
Imaging
No	—	—
Yes	0.38	0.07, 1.60	0.2
Centre Volume	1.01	0.96, 1.06	0.6
Age at Time of Injury	1.00	0.95, 1.06	>0.9
Urethral Injury
No	—	—
Yes	0.88	0.16, 3.60	0.9
Transfer for Surgical Repair
No	—	—
Yes	0.73	0.13, 3.47	0.7
Site of Injury
Crura	—	—
Distal	0.56	0.01, 123	0.8
Mid	0.99	0.03, 194	>0.9
Multi-level	0.76	0.00, 223	>0.9
Proximal	1.02	0.03, 197	>0.9
Side of Injury
Bilateral	—	—
Unilateral	0.34	0.08, 1.48	0.15
Mechanism of Injury
Sexual intercourse trauma	—	—
Other	0.48	0.00, 7.84	0.7
Non-sexual blunt force trauma	0.18	0.00, 1.91	0.2
Masturbation injury	0.21	0.00, 2.65	0.3
Abbreviation: CI = Confidence Interval
Number of patients in analysis: 119

7.3.3 ROC

waisting_shortening_lr_predict <-
  predict(waisting_shortening_lr, waisting_shortening_for_log_reg_18, type = "response")

waisting_shortening_outcome <- as.factor(ifelse(waisting_shortening_for_log_reg_18$waisting_shortening == "Yes",
                                  1,
                                  0))

waisting_shortening_lr_roc <-
  pROC::roc(waisting_shortening_outcome, waisting_shortening_lr_predict)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

auc(waisting_shortening_lr_roc)

## Area under the curve: 0.7673

ci.auc(waisting_shortening_lr_roc)

## 95% CI: 0.6452-0.8894 (DeLong)

waisting_shortening_lr_roc_data <- cbind(sensitivities = waisting_shortening_lr_roc$sensitivities,
                                specificities = waisting_shortening_lr_roc$specificities) %>% as_tibble()

ggplot(as_tibble(waisting_shortening_lr_roc_data), aes(x = 1 - specificities, y = sensitivities)) +
  geom_line() +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed")

7.3.4 Calibration plot

waisting_shortening_lr_predict2 <- as_tibble(waisting_shortening_lr_predict) %>% cbind(waisting_shortening_for_log_reg_18) %>% subset(select = c(
  value,
  waisting_shortening,
  did_the_patient_get_to_theatre_within_18_hours_of_presentation
)) %>% mutate(
  pred = value,
  .keep = "unused"
)

waisting_shortening_lr_predict2$waisting_shortening <- ifelse(waisting_shortening_lr_predict2$waisting_shortening == "Yes", 1, 0)

waisting_shortening_lr_predict2$did_the_patient_get_to_theatre_within_18_hours_of_presentation <- as.integer(waisting_shortening_lr_predict2$did_the_patient_get_to_theatre_within_18_hours_of_presentation)

## Warning: NAs introduced by coercion

predtools::calibration_plot(
  data = waisting_shortening_lr_predict2,
  obs = "waisting_shortening",
  pred = "pred",
  y_lim = c(0, 1),
  x_lim = c(0, 1),
  title = "Calibration plot for waisting_shortening Model"
)

## $calibration_plot

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

7.3.5 Brier Score

brier <- function(model, data, outcome) {
  probs <- predict(model, type = "response")
  obs <- data[[outcome]]
  mean((obs - probs)^2, na.rm = TRUE)
}

brier_score <- brier(waisting_shortening_lr, waisting_shortening_for_log_reg_182, "waisting_shortening") 

cbind(
  "Score" = "Brier Score",
  "Brier Score" = round(brier_score, digits = 3)
) %>% as_tibble() %>% gt() %>% gt_theme_espn()

Score	Brier Score
Brier Score	0.1

7.3.6 Nagelkerke’s R squared

pseudoR2_logistf(waisting_shortening_lr, waisting_shortening_for_log_reg_182)

Score	Nagelkerke's R^2
Nagelkerke's R^2	0.157

7.4 Abnormal Curvature

7.4.1 Post-hoc Power calculation

n1 <- nrow(abnormal_curvature_for_log_reg_18 %>% subset(did_the_patient_get_to_theatre_within_18_hours_of_presentation == "No"))
p1 <- nrow(abnormal_curvature_for_log_reg_18 %>% subset(did_the_patient_get_to_theatre_within_18_hours_of_presentation == "No" & abnormal_curvature == "Yes")) / n1

n2 <- nrow(abnormal_curvature_for_log_reg_18 %>% subset(did_the_patient_get_to_theatre_within_18_hours_of_presentation == "Yes"))
alpha <- 0.05

p2 <- nrow(abnormal_curvature_for_log_reg_18 %>% subset(did_the_patient_get_to_theatre_within_18_hours_of_presentation == "Yes" & abnormal_curvature == "Yes")) / n2

power <- calculate_power_table(p1, p2, n1, n2, alpha, title = "Power Calculation Results for Post-Operative abnormal_curvature")
power

Power Calculation Results for Post-Operative abnormal_curvature
Based on Two Proportion Test
Parameter	Value
Group 1 Proportion (p1)	0.18
Group 2 Proportion (p2)	0.14
Group 1 Size (n1)	112.00
Group 2 Size (n2)	133.00
Significance Level (alpha)	0.05
Calculated Power (%)	0.30

7.4.2 Logistic Regression (Firth Correction)

Multi-level removed as model failed to separate

abnormal_curvature_for_log_reg_18 <- abnormal_curvature_for_log_reg_18 %>% 
  left_join(knife_to_skin1,
            by = c("submission_id" = "submission_id")) %>%
  mutate(abnormal_curvature = abnormal_curvature %>% factor(levels = c("Yes",
                                           "No")))

abnormal_curvature_for_log_reg_182 <- abnormal_curvature_for_log_reg_18 %>%
  mutate(presentation_to_knife_to_skin = ifelse(
    presentation_to_knife_to_skin < 0 | presentation_to_knife_to_skin > 5000,
    NA,
    presentation_to_knife_to_skin / (60*6)
  ),
    abnormal_curvature = case_when(abnormal_curvature == "Yes" ~ 1,
                     abnormal_curvature == "No" ~ 0,
                     TRUE ~ NA_real_),
  .keep = "unused"
  ) %>%
  select(
    abnormal_curvature,
    imaging,
    volume_n,
    what_age_was_the_patient_at_the_time_of_the_injury,
    did_the_patient_get_to_theatre_within_18_hours_of_presentation,
    what_age_was_the_patient_at_the_time_of_the_injury,
    was_the_patient_transferred_from_another_hospital_for_surgical_repair,
    was_a_urethral_injury_noted_at_surgery,
    laterality,
    distal,
    proximal,
    mid,
    crura,
    multi_level,
    mechanism_of_injury
  ) %>% drop_na()

abnormal_curvature_lr <- logistf(
  abnormal_curvature ~ did_the_patient_get_to_theatre_within_18_hours_of_presentation + imaging + volume_n + 
    what_age_was_the_patient_at_the_time_of_the_injury + was_a_urethral_injury_noted_at_surgery +
    was_the_patient_transferred_from_another_hospital_for_surgical_repair + distal + proximal + mid +
    crura + laterality + mechanism_of_injury,
  data = abnormal_curvature_for_log_reg_182
)



abnormal_curvature_lr %>%
  tbl_regression(
    exponentiate = TRUE,
    label = list(
      did_the_patient_get_to_theatre_within_18_hours_of_presentation = "Theatre within 18 hours of Presentation",
      was_a_urethral_injury_noted_at_surgery = "Urethral Injury",
      was_the_patient_transferred_from_another_hospital_for_surgical_repair = "Transfer for Surgical Repair",
      locality = "Site of Injury",
      distal = "Distal Corpus Location",
    proximal = "Proximal Corpus Location",
    mid = "Mid Corpus Location",
    crura = "Crural location",
      mechanism_of_injury = "Mechanism of Injury",
      what_age_was_the_patient_at_the_time_of_the_injury =
        "Age at Time of Injury",
    imaging = "Imaging",
    volume_n = "Centre Volume"
    )
  ) %>%
  bold_labels() %>%
  as_gt() %>%  
  tab_source_note(
    md(glue(
      "Number of patients in analysis: {nrow(abnormal_curvature_for_log_reg_182)}"
    ))
  )

## ! `broom::tidy()` failed to tidy the model.

## ✔ `tidy_parameters()` used instead.

## ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.

Characteristic	exp(Beta)	95% CI	p-value
Theatre within 18 hours of Presentation
No	—	—
Yes	0.65	0.25, 1.63	0.4
Imaging
No	—	—
Yes	0.37	0.12, 1.10	0.073
Centre Volume	1.03	0.99, 1.06	0.2
Age at Time of Injury	0.97	0.93, 1.01	0.15
Urethral Injury
No	—	—
Yes	1.10	0.36, 3.11	0.9
Transfer for Surgical Repair
No	—	—
Yes	1.77	0.59, 5.34	0.3
Distal Corpus Location
No	—	—
Yes	0.42	0.00, 5.93	0.6
Proximal Corpus Location
No	—	—
Yes	0.40	0.00, 6.38	0.6
Mid Corpus Location
No	—	—
Yes	0.24	0.00, 3.43	0.3
Crural location
No	—	—
Yes	0.09	0.00, 7.45	0.3
laterality
Bilateral	—	—
Unilateral	0.56	0.16, 1.93	0.4
Mechanism of Injury
Sexual intercourse trauma	—	—
Other	0.28	0.00, 3.29	0.4
Non-sexual blunt force trauma	0.44	0.09, 1.76	0.3
Masturbation injury	0.12	0.00, 1.09	0.062
Abbreviation: CI = Confidence Interval
Number of patients in analysis: 128

7.4.3 ROC

abnormal_curvature_lr_predict <-
  predict(abnormal_curvature_lr, abnormal_curvature_for_log_reg_18, type = "response")

abnormal_curvature_outcome <- as.factor(ifelse(abnormal_curvature_for_log_reg_18$abnormal_curvature == "Yes",
                                  1,
                                  0))

abnormal_curvature_lr_roc <-
  pROC::roc(abnormal_curvature_outcome, abnormal_curvature_lr_predict)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

auc(abnormal_curvature_lr_roc)

## Area under the curve: 0.7663

ci.auc(abnormal_curvature_lr_roc)

## 95% CI: 0.6832-0.8493 (DeLong)

abnormal_curvature_lr_roc_data <- cbind(sensitivities = abnormal_curvature_lr_roc$sensitivities,
                                specificities = abnormal_curvature_lr_roc$specificities) %>% as_tibble()

ggplot(as_tibble(abnormal_curvature_lr_roc_data), aes(x = 1 - specificities, y = sensitivities)) +
  geom_line() +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed")

7.4.4 Calibration plot

abnormal_curvature_lr_predict2 <- as_tibble(abnormal_curvature_lr_predict) %>% cbind(abnormal_curvature_for_log_reg_18) %>% subset(select = c(
  value,
  abnormal_curvature,
  did_the_patient_get_to_theatre_within_18_hours_of_presentation
)) %>% mutate(
  pred = value,
  .keep = "unused"
)

abnormal_curvature_lr_predict2$abnormal_curvature <- ifelse(abnormal_curvature_lr_predict2$abnormal_curvature == "Yes", 1, 0)

abnormal_curvature_lr_predict2$did_the_patient_get_to_theatre_within_18_hours_of_presentation <- as.integer(abnormal_curvature_lr_predict2$did_the_patient_get_to_theatre_within_18_hours_of_presentation)

## Warning: NAs introduced by coercion

predtools::calibration_plot(
  data = abnormal_curvature_lr_predict2,
  obs = "abnormal_curvature",
  pred = "pred",
  y_lim = c(0, 1),
  x_lim = c(0, 1),
  title = "Calibration plot for abnormal_curvature Model"
)

## $calibration_plot

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

7.4.5 Brier Score

brier <- function(model, data, outcome) {
  probs <- predict(model, type = "response")
  obs <- data[[outcome]]
  mean((obs - probs)^2, na.rm = TRUE)
}

brier_score <- brier(abnormal_curvature_lr, abnormal_curvature_for_log_reg_182, "abnormal_curvature") 

cbind(
  "Score" = "Brier Score",
  "Brier Score" = round(brier_score, digits = 3)
) %>% as_tibble() %>% gt() %>% gt_theme_espn()

Score	Brier Score
Brier Score	0.171

7.4.6 Nagelkerke’s R squared

pseudoR2_logistf(abnormal_curvature_lr, abnormal_curvature_for_log_reg_182)

Score	Nagelkerke's R^2
Nagelkerke's R^2	0.229

8 Logistic Regression - Time from Presentation to Theatre with 36 hour cut-off

8.1 Pain

8.1.1 Post-hoc Power calculation

n1 <- nrow(pain_for_log_reg_36 %>% subset(did_the_patient_get_to_theatre_within_36_hours_of_presentation == "No"))
p1 <- nrow(pain_for_log_reg_36 %>% subset(did_the_patient_get_to_theatre_within_36_hours_of_presentation == "No" & pain == "Yes")) / n1

n2 <- nrow(pain_for_log_reg_36 %>% subset(did_the_patient_get_to_theatre_within_36_hours_of_presentation == "Yes"))
alpha <- 0.05

p2 <- nrow(pain_for_log_reg_36 %>% subset(did_the_patient_get_to_theatre_within_36_hours_of_presentation == "Yes" & pain == "Yes")) / n2

power <- calculate_power_table(p1, p2, n1, n2, alpha, title = "Power Calculation Results for Post-Operative Pain")
power

Power Calculation Results for Post-Operative Pain
Based on Two Proportion Test
Parameter	Value
Group 1 Proportion (p1)	0.07
Group 2 Proportion (p2)	0.12
Group 1 Size (n1)	45.00
Group 2 Size (n2)	200.00
Significance Level (alpha)	0.05
Calculated Power (%)	26.10

8.1.2 Logistic Regression (Firth Correction)

knife_to_skin1 <- knife_to_skin %>% 
  subset(select = c(
    submission_id,
    time_from_injury_to_th
  ))

pain_for_log_reg_36 <- pain_for_log_reg_36 %>% 
  left_join(knife_to_skin1,
            by = c("submission_id" = "submission_id")) %>%
  mutate(pain = pain %>% factor(levels = c("Yes",
                                           "No")))

pain_for_log_reg_362 <- pain_for_log_reg_36 %>%
  mutate(presentation_to_knife_to_skin = ifelse(
    presentation_to_knife_to_skin < 0 | presentation_to_knife_to_skin > 7920,
    NA,
    presentation_to_knife_to_skin / (60*6)
  ),
    pain = case_when(pain == "Yes" ~ 1,
                     pain == "No" ~ 0,
                     TRUE ~ NA_real_),
  .keep = "unused"
  ) %>%
  select(
    pain,
    imaging,
    volume_n,
    what_age_was_the_patient_at_the_time_of_the_injury,
    did_the_patient_get_to_theatre_within_36_hours_of_presentation,
    what_age_was_the_patient_at_the_time_of_the_injury,
    was_the_patient_transferred_from_another_hospital_for_surgical_repair,
    was_a_urethral_injury_noted_at_surgery,
    laterality,
    distal,
    proximal,
    mid,
    crura,
    multi_level,
    mechanism_of_injury
  ) %>% drop_na()

pain_lr <- logistf(
  pain ~ did_the_patient_get_to_theatre_within_36_hours_of_presentation + imaging + volume_n +
    what_age_was_the_patient_at_the_time_of_the_injury + was_a_urethral_injury_noted_at_surgery +
    was_the_patient_transferred_from_another_hospital_for_surgical_repair + distal + proximal + mid +
    crura + multi_level + laterality + mechanism_of_injury,
  data = pain_for_log_reg_362
)

pain_lr %>%
  tbl_regression(
    exponentiate = TRUE,
    label = list(
      did_the_patient_get_to_theatre_within_36_hours_of_presentation = "Theatre within 36 hours of Presentation",
      was_a_urethral_injury_noted_at_surgery = "Urethral Injury",
      was_the_patient_transferred_from_another_hospital_for_surgical_repair = "Transfer for Surgical Repair",
      laterality = "Side of Injury",
      distal = "Distal Corpus Location",
    proximal = "Proximal Corpus Location",
    mid = "Mid Corpus Location",
    crura = "Crural location",
    multi_level = "Multi-level Injury",
      mechanism_of_injury = "Mechanism of Injury",
      what_age_was_the_patient_at_the_time_of_the_injury =
        "Age at Time of Injury",
    imaging = "Imaging",
    volume_n = "Centre Volume"
    )
  ) %>%
  bold_labels() %>%
  as_gt() %>%  
  tab_source_note(
    md(glue(
      "Number of patients in analysis: {nrow(pain_for_log_reg_362)}"
    ))
  )

## ! `broom::tidy()` failed to tidy the model.

## ✔ `tidy_parameters()` used instead.

## ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.

Characteristic	exp(Beta)	95% CI	p-value
Theatre within 36 hours of Presentation
No	—	—
Yes	2.36	0.66, 11.3	0.2
Imaging
No	—	—
Yes	1.00	0.36, 2.88	>0.9
Centre Volume	0.94	0.89, 0.98	0.006
Age at Time of Injury	1.04	0.99, 1.09	0.13
Urethral Injury
No	—	—
Yes	1.22	0.34, 3.94	0.7
Transfer for Surgical Repair
No	—	—
Yes	4.06	1.26, 14.0	0.019
Distal Corpus Location
No	—	—
Yes	0.01	0.00, 28.3	0.3
Proximal Corpus Location
No	—	—
Yes	0.01	0.00, 27.1	0.3
Mid Corpus Location
No	—	—
Yes	0.01	0.00, 15.6	0.2
Crural location
No	—	—
Yes	0.15	0.00, 75.7	0.5
Multi-level Injury
No	—	—
Yes	508	0.07, 8,497,815	0.2
Side of Injury
Bilateral	—	—
Unilateral	0.64	0.16, 2.91	0.5
Mechanism of Injury
Sexual intercourse trauma	—	—
Other	0.69	0.00, 9.64	0.8
Non-sexual blunt force trauma	3.47	0.71, 15.4	0.12
Masturbation injury	4.62	0.67, 32.7	0.11
Abbreviation: CI = Confidence Interval
Number of patients in analysis: 141

# Add back in to get 24hr cut-off: presentation_to_knife_to_skin,
# did_the_patient_get_to_theatre_within_36_hours_of_presentation
# time_from_injury_to_th

8.1.3 ROC

pain_lr_predict <-
  predict(pain_lr, pain_for_log_reg_36, type = "response")

pain_outcome <- as.factor(ifelse(pain_for_log_reg_36$pain == "Yes",
                                  1,
                                  0))

pain_lr_roc <-
  pROC::roc(pain_outcome, pain_lr_predict)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

auc(pain_lr_roc)

## Area under the curve: 0.7904

ci.auc(pain_lr_roc)

## 95% CI: 0.7005-0.8804 (DeLong)

pain_lr_roc_data <- cbind(sensitivities = pain_lr_roc$sensitivities,
                                specificities = pain_lr_roc$specificities) %>% as_tibble()

ggplot(as_tibble(pain_lr_roc_data), aes(x = 1 - specificities, y = sensitivities)) +
  geom_line() +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed")

8.1.4 Calibration plot

pain_lr_predict2 <- as_tibble(pain_lr_predict) %>% cbind(pain_for_log_reg_36) %>% subset(select = c(
  value,
  pain,
  did_the_patient_get_to_theatre_within_36_hours_of_presentation
)) %>% mutate(
  pred = value,
  .keep = "unused"
)

pain_lr_predict2$pain <- ifelse(pain_lr_predict2$pain == "Yes", 1, 0)

pain_lr_predict2$did_the_patient_get_to_theatre_within_36_hours_of_presentation <- as.integer(pain_lr_predict2$did_the_patient_get_to_theatre_within_36_hours_of_presentation)

## Warning: NAs introduced by coercion

predtools::calibration_plot(
  data = pain_lr_predict2,
  obs = "pain",
  pred = "pred",
  y_lim = c(0, 1),
  x_lim = c(0, 1),
  title = "Calibration plot for Pain Model"
)

## $calibration_plot

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

8.1.5 Brier Score

brier <- function(model, data, outcome) {
  probs <- predict(model, type = "response")
  obs <- data[[outcome]]
  mean((obs - probs)^2, na.rm = TRUE)
}

brier_score <- brier(pain_lr, pain_for_log_reg_362, "pain") 

cbind(
  "Score" = "Brier Score",
  "Brier Score" = round(brier_score, digits = 3)
) %>% as_tibble() %>% gt() %>% gt_theme_espn()

Score	Brier Score
Brier Score	0.117

8.1.6 Nagelkerke’s R squared

pseudoR2_logistf <- function(model, data) {
  # Extract log-likelihoods
  ll_full <- model$loglik["full"]
  ll_null <- model$loglik["null"]
  
  # Sample size
  n <- nrow(data)
  
  # Cox & Snell R²
  R2_CS <- 1 - exp((2 / n) * (ll_null - ll_full))
  
  # Nagelkerke R²
  R2_N <- R2_CS / (1 - exp((2 / n) * ll_null))

  R2_N <- R2_N %>% as_tibble()
  
  x <- cbind(
  "Score" = "Nagelkerke's R^2",
  "Nagelkerke's R^2" = round(R2_N$value, digits = 3)
) %>% as_tibble() %>% gt() %>% gt_theme_espn()
  
  # Return as a named list
  return(
    x
  )
}

pseudoR2_logistf(pain_lr, pain_for_log_reg_362)

Score	Nagelkerke's R^2
Nagelkerke's R^2	0.263

8.2 ED

8.2.1 Post-hoc Power calculation

n1 <- nrow(ed_for_log_reg_36 %>% subset(did_the_patient_get_to_theatre_within_36_hours_of_presentation == "No"))
p1 <- nrow(ed_for_log_reg_36 %>% subset(did_the_patient_get_to_theatre_within_36_hours_of_presentation == "No" & ed == "Yes")) / n1

n2 <- nrow(ed_for_log_reg_36 %>% subset(did_the_patient_get_to_theatre_within_36_hours_of_presentation == "Yes"))
alpha <- 0.05

p2 <- nrow(ed_for_log_reg_36 %>% subset(did_the_patient_get_to_theatre_within_36_hours_of_presentation == "Yes" & ed == "Yes")) / n2

power <- calculate_power_table(p1, p2, n1, n2, alpha, title = "Power Calculation Results for Post-Operative Erectile Dysfunction")
power

Power Calculation Results for Post-Operative Erectile Dysfunction
Based on Two Proportion Test
Parameter	Value
Group 1 Proportion (p1)	0.27
Group 2 Proportion (p2)	0.14
Group 1 Size (n1)	45.00
Group 2 Size (n2)	200.00
Significance Level (alpha)	0.05
Calculated Power (%)	0.00

8.2.2 Logistic Regression (Firth Correction)

ed_for_log_reg_36 <- ed_for_log_reg_36 %>%
  left_join(knife_to_skin1, by = c("submission_id" = "submission_id")) %>%
  mutate(ed = ed %>% factor(levels = c("No", "Yes")))

ed_for_log_reg_362 <- ed_for_log_reg_36 %>%
  drop_na(
    ed,
    imaging,
    volume_n,
    did_the_patient_get_to_theatre_within_36_hours_of_presentation,
    what_age_was_the_patient_at_the_time_of_the_injury,
    presentation_to_knife_to_skin,
    what_age_was_the_patient_at_the_time_of_the_injury,
    was_the_patient_transferred_from_another_hospital_for_surgical_repair,
    was_a_urethral_injury_noted_at_surgery,
    laterality,
    distal,
    proximal,
    mid,
    crura,
    multi_level,
    mechanism_of_injury
  ) %>%
  mutate(
    ed = case_when(ed == "Yes" ~ 1,
                   ed == "No" ~ 0,
                   TRUE ~ NA_real_),
    presentation_to_knife_to_skin = ifelse(
    presentation_to_knife_to_skin < 0 | presentation_to_knife_to_skin > 5000,
    NA,
    presentation_to_knife_to_skin / (60*6)
  ),
    .keep = "unused"
  )

ed_lr <- logistf(
  ed ~ did_the_patient_get_to_theatre_within_36_hours_of_presentation + imaging + volume_n +
    what_age_was_the_patient_at_the_time_of_the_injury + was_a_urethral_injury_noted_at_surgery +
    was_the_patient_transferred_from_another_hospital_for_surgical_repair + distal + proximal + mid +
    crura + multi_level + laterality + mechanism_of_injury,
  data = ed_for_log_reg_362
)

ed_lr %>%
  tbl_regression(
    exponentiate = TRUE,
    label = list(
      did_the_patient_get_to_theatre_within_36_hours_of_presentation = "Theatre within 36 hours of Presentation",
      was_a_urethral_injury_noted_at_surgery = "Urethral Injury",
      was_the_patient_transferred_from_another_hospital_for_surgical_repair = "Transfer for Surgical Repair",
      distal = "Distal Corpus Location",
    proximal = "Proximal Corpus Location",
    mid = "Mid Corpus Location",
    crura = "Crural location",
    multi_level = "Multi-level Injury",
      laterality = "Side of Injury",
      mechanism_of_injury = "Mechanism of Injury",
      what_age_was_the_patient_at_the_time_of_the_injury =
        "Age at Time of Injury",
    imaging = "Imaging",
    volume_n = "Centre Volume"
    )
  ) %>%
  bold_labels() %>%
  as_gt() %>%  
  tab_source_note(
    md(glue(
      "Number of patients in analysis: {nrow(ed_for_log_reg_362)}"
    ))
  )

## ! `broom::tidy()` failed to tidy the model.

## ✔ `tidy_parameters()` used instead.

## ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.

Characteristic	exp(Beta)	95% CI	p-value
Theatre within 36 hours of Presentation
No	—	—
Yes	0.37	0.13, 1.03	0.057
Imaging
No	—	—
Yes	1.11	0.39, 3.24	0.8
Centre Volume	0.98	0.94, 1.02	0.3
Age at Time of Injury	1.02	0.98, 1.06	0.4
Urethral Injury
No	—	—
Yes	3.03	1.00, 9.63	0.051
Transfer for Surgical Repair
No	—	—
Yes	1.88	0.60, 6.25	0.3
Distal Corpus Location
No	—	—
Yes	40.0	0.06, 108,356	0.3
Proximal Corpus Location
No	—	—
Yes	27.3	0.05, 58,756	0.3
Mid Corpus Location
No	—	—
Yes	11.7	0.02, 25,461	0.5
Crural location
No	—	—
Yes	12.9	0.06, 11,323	0.4
Multi-level Injury
No	—	—
Yes	0.07	0.00, 167	0.5
Side of Injury
Bilateral	—	—
Unilateral	7.04	1.45, 53.2	0.013
Mechanism of Injury
Sexual intercourse trauma	—	—
Other	0.28	0.00, 4.42	0.4
Non-sexual blunt force trauma	0.29	0.03, 1.43	0.14
Masturbation injury	0.12	0.00, 1.26	0.085
Abbreviation: CI = Confidence Interval
Number of patients in analysis: 136

# did_the_patient_get_to_theatre_within_36_hours_of_presentation,

8.2.3 ROC

ed_lr_predict <-
  predict(ed_lr, ed_for_log_reg_36, type = "response")

ed_outcome <- as.factor(ifelse(ed_for_log_reg_36$ed == "Yes",
                                  1,
                                  0))

ed_lr_roc <-
  pROC::roc(ed_outcome, ed_lr_predict)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

auc(ed_lr_roc)

## Area under the curve: 0.8143

ci.auc(ed_lr_roc)

## 95% CI: 0.7345-0.8941 (DeLong)

ed_lr_roc_data <- cbind(sensitivities = ed_lr_roc$sensitivities,
                                specificities = ed_lr_roc$specificities) %>% as_tibble()

ggplot(as_tibble(ed_lr_roc_data), aes(x = 1 - specificities, y = sensitivities)) +
  geom_line() +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed")

8.2.4 Calibration plot

ed_lr_predict2 <- as_tibble(ed_lr_predict) %>% cbind(ed_for_log_reg_36) %>% subset(select = c(
  value,
  ed,
  did_the_patient_get_to_theatre_within_36_hours_of_presentation
)) %>% mutate(
  pred = value,
  .keep = "unused"
)

ed_lr_predict2$ed <- ifelse(ed_lr_predict2$ed == "Yes", 1, 0)

ed_lr_predict2$did_the_patient_get_to_theatre_within_36_hours_of_presentation <- as.integer(ed_lr_predict2$did_the_patient_get_to_theatre_within_36_hours_of_presentation)

## Warning: NAs introduced by coercion

predtools::calibration_plot(
  data = ed_lr_predict2,
  obs = "ed",
  pred = "pred",
  y_lim = c(0, 1),
  x_lim = c(0, 1),
  title = "Calibration plot for ed Model"
)

## $calibration_plot

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

8.2.5 Brier Score

brier_score <- brier(ed_lr, ed_for_log_reg_362, "ed")

cbind(
  "Score" = "Brier Score",
  "Brier Score" = round(brier_score, digits = 3)
) %>% as_tibble() %>% gt() %>% gt_theme_espn()

Score	Brier Score
Brier Score	0.155

8.2.6 Nagelkerke’s R squared

pseudoR2_logistf(ed_lr, ed_for_log_reg_362)

Score	Nagelkerke's R^2
Nagelkerke's R^2	0.265

8.3 Waisting / Shortening

8.3.1 Post-hoc Power calculation

n1 <- nrow(waisting_shortening_for_log_reg_36 %>% subset(did_the_patient_get_to_theatre_within_36_hours_of_presentation == "No"))
p1 <- nrow(waisting_shortening_for_log_reg_36 %>% subset(did_the_patient_get_to_theatre_within_36_hours_of_presentation == "No" & waisting_shortening == "Yes")) / n1

n2 <- nrow(waisting_shortening_for_log_reg_36 %>% subset(did_the_patient_get_to_theatre_within_36_hours_of_presentation == "Yes"))
alpha <- 0.05

p2 <- nrow(waisting_shortening_for_log_reg_36 %>% subset(did_the_patient_get_to_theatre_within_36_hours_of_presentation == "Yes" & waisting_shortening == "Yes")) / n2

power <- calculate_power_table(p1, p2, n1, n2, alpha, title = "Power Calculation Results for Post-Operative waisting_shortening")
power

Power Calculation Results for Post-Operative waisting_shortening
Based on Two Proportion Test
Parameter	Value
Group 1 Proportion (p1)	0.07
Group 2 Proportion (p2)	0.07
Group 1 Size (n1)	45.00
Group 2 Size (n2)	200.00
Significance Level (alpha)	0.05
Calculated Power (%)	3.00

8.3.2 Logistic Regression (Firth Correction)

waisting_shortening_for_log_reg_36 <- waisting_shortening_for_log_reg_36 %>% 
  left_join(knife_to_skin1,
            by = c("submission_id" = "submission_id")) %>%
  mutate(waisting_shortening = waisting_shortening %>% factor(levels = c("Yes",
                                           "No")))

waisting_shortening_for_log_reg_362 <- waisting_shortening_for_log_reg_36 %>%
  mutate(presentation_to_knife_to_skin = ifelse(
    presentation_to_knife_to_skin < 0 | presentation_to_knife_to_skin > 5000,
    NA,
    presentation_to_knife_to_skin / (60*6)
  ),
    waisting_shortening = case_when(waisting_shortening == "Yes" ~ 1,
                     waisting_shortening == "No" ~ 0,
                     TRUE ~ NA_real_),
  .keep = "unused"
  ) %>%
  select(
    waisting_shortening,
    imaging,
    volume_n,
    what_age_was_the_patient_at_the_time_of_the_injury,
    did_the_patient_get_to_theatre_within_36_hours_of_presentation,
    what_age_was_the_patient_at_the_time_of_the_injury,
    was_the_patient_transferred_from_another_hospital_for_surgical_repair,
    was_a_urethral_injury_noted_at_surgery,
    laterality,
    locality,
    mechanism_of_injury
  ) %>% drop_na()

waisting_shortening_lr <- logistf(
  waisting_shortening ~ did_the_patient_get_to_theatre_within_36_hours_of_presentation + imaging + volume_n +
    what_age_was_the_patient_at_the_time_of_the_injury + was_a_urethral_injury_noted_at_surgery +
    was_the_patient_transferred_from_another_hospital_for_surgical_repair + locality + laterality + mechanism_of_injury,
  data = waisting_shortening_for_log_reg_362
)

waisting_shortening_lr %>%
  tbl_regression(
    exponentiate = TRUE,
    label = list(
      did_the_patient_get_to_theatre_within_36_hours_of_presentation = "Theatre within 36 hours of Presentation",
      was_a_urethral_injury_noted_at_surgery = "Urethral Injury",
      was_the_patient_transferred_from_another_hospital_for_surgical_repair = "Transfer for Surgical Repair",
      locality = "Site of Injury",
      laterality = "Side of Injury",
      mechanism_of_injury = "Mechanism of Injury",
      what_age_was_the_patient_at_the_time_of_the_injury =
        "Age at Time of Injury",
      imaging = "Imaging",
      volume_n = "Centre Volume"
    )
  ) %>%
  bold_labels() %>%
  as_gt() %>%  
  tab_source_note(
    md(glue(
      "Number of patients in analysis: {nrow(waisting_shortening_for_log_reg_362)}"
    ))
  )

## ! `broom::tidy()` failed to tidy the model.

## ✔ `tidy_parameters()` used instead.

## ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.

Characteristic	exp(Beta)	95% CI	p-value
Theatre within 36 hours of Presentation
No	—	—
Yes	1.12	0.28, 5.30	0.9
Imaging
No	—	—
Yes	0.59	0.14, 2.37	0.5
Centre Volume	1.01	0.96, 1.06	0.6
Age at Time of Injury	1.01	0.95, 1.07	0.7
Urethral Injury
No	—	—
Yes	0.96	0.18, 3.78	>0.9
Transfer for Surgical Repair
No	—	—
Yes	0.61	0.11, 2.82	0.5
Site of Injury
Crura	—	—
Distal	0.43	0.01, 94.3	0.7
Mid	0.76	0.02, 145	0.9
Multi-level	0.45	0.00, 123	0.7
Proximal	0.67	0.02, 127	0.8
Side of Injury
Bilateral	—	—
Unilateral	0.30	0.07, 1.25	0.10
Mechanism of Injury
Sexual intercourse trauma	—	—
Other	0.48	0.00, 6.88	0.6
Non-sexual blunt force trauma	0.24	0.00, 2.15	0.2
Masturbation injury	0.32	0.00, 3.53	0.4
Abbreviation: CI = Confidence Interval
Number of patients in analysis: 119

8.3.3 ROC

waisting_shortening_lr_predict <-
  predict(waisting_shortening_lr, waisting_shortening_for_log_reg_36, type = "response")

waisting_shortening_outcome <- as.factor(ifelse(waisting_shortening_for_log_reg_36$waisting_shortening == "Yes",
                                  1,
                                  0))

waisting_shortening_lr_roc <-
  pROC::roc(waisting_shortening_outcome, waisting_shortening_lr_predict)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

auc(waisting_shortening_lr_roc)

## Area under the curve: 0.7667

ci.auc(waisting_shortening_lr_roc)

## 95% CI: 0.6419-0.8915 (DeLong)

waisting_shortening_lr_roc_data <- cbind(sensitivities = waisting_shortening_lr_roc$sensitivities,
                                specificities = waisting_shortening_lr_roc$specificities) %>% as_tibble()

ggplot(as_tibble(waisting_shortening_lr_roc_data), aes(x = 1 - specificities, y = sensitivities)) +
  geom_line() +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed")

8.3.4 Calibration plot

waisting_shortening_lr_predict2 <- as_tibble(waisting_shortening_lr_predict) %>% cbind(waisting_shortening_for_log_reg_36) %>% subset(select = c(
  value,
  waisting_shortening,
  did_the_patient_get_to_theatre_within_36_hours_of_presentation
)) %>% mutate(
  pred = value,
  .keep = "unused"
)

waisting_shortening_lr_predict2$waisting_shortening <- ifelse(waisting_shortening_lr_predict2$waisting_shortening == "Yes", 1, 0)

waisting_shortening_lr_predict2$did_the_patient_get_to_theatre_within_36_hours_of_presentation <- as.integer(waisting_shortening_lr_predict2$did_the_patient_get_to_theatre_within_36_hours_of_presentation)

## Warning: NAs introduced by coercion

predtools::calibration_plot(
  data = waisting_shortening_lr_predict2,
  obs = "waisting_shortening",
  pred = "pred",
  y_lim = c(0, 1),
  x_lim = c(0, 1),
  title = "Calibration plot for waisting_shortening Model"
)

## $calibration_plot

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

8.3.5 Brier Score

brier <- function(model, data, outcome) {
  probs <- predict(model, type = "response")
  obs <- data[[outcome]]
  mean((obs - probs)^2, na.rm = TRUE)
}

brier_score <- brier(waisting_shortening_lr, waisting_shortening_for_log_reg_362, "waisting_shortening") 

cbind(
  "Score" = "Brier Score",
  "Brier Score" = round(brier_score, digits = 3)
) %>% as_tibble() %>% gt() %>% gt_theme_espn()

Score	Brier Score
Brier Score	0.102

8.3.6 Nagelkerke’s R squared

pseudoR2_logistf(waisting_shortening_lr, waisting_shortening_for_log_reg_362)

Score	Nagelkerke's R^2
Nagelkerke's R^2	0.128

8.4 Abnormal Curvature

8.4.1 Post-hoc Power calculation

n1 <- nrow(abnormal_curvature_for_log_reg_36 %>% subset(did_the_patient_get_to_theatre_within_36_hours_of_presentation == "No"))
p1 <- nrow(abnormal_curvature_for_log_reg_36 %>% subset(did_the_patient_get_to_theatre_within_36_hours_of_presentation == "No" & abnormal_curvature == "Yes")) / n1

n2 <- nrow(abnormal_curvature_for_log_reg_36 %>% subset(did_the_patient_get_to_theatre_within_36_hours_of_presentation == "Yes"))
alpha <- 0.05

p2 <- nrow(abnormal_curvature_for_log_reg_36 %>% subset(did_the_patient_get_to_theatre_within_36_hours_of_presentation == "Yes" & abnormal_curvature == "Yes")) / n2

power <- calculate_power_table(p1, p2, n1, n2, alpha, title = "Power Calculation Results for Post-Operative abnormal_curvature")
power

Power Calculation Results for Post-Operative abnormal_curvature
Based on Two Proportion Test
Parameter	Value
Group 1 Proportion (p1)	0.20
Group 2 Proportion (p2)	0.15
Group 1 Size (n1)	45.00
Group 2 Size (n2)	200.00
Significance Level (alpha)	0.05
Calculated Power (%)	0.30

8.4.2 Logistic Regression (Firth Correction)

Multi-level removed as model failed to separate

abnormal_curvature_for_log_reg_36 <- abnormal_curvature_for_log_reg_36 %>% 
  left_join(knife_to_skin1,
            by = c("submission_id" = "submission_id")) %>%
  mutate(abnormal_curvature = abnormal_curvature %>% factor(levels = c("Yes",
                                           "No")))

abnormal_curvature_for_log_reg_362 <- abnormal_curvature_for_log_reg_36 %>%
  mutate(presentation_to_knife_to_skin = ifelse(
    presentation_to_knife_to_skin < 0 | presentation_to_knife_to_skin > 5000,
    NA,
    presentation_to_knife_to_skin / (60*6)
  ),
    abnormal_curvature = case_when(abnormal_curvature == "Yes" ~ 1,
                     abnormal_curvature == "No" ~ 0,
                     TRUE ~ NA_real_),
  .keep = "unused"
  ) %>%
  select(
    abnormal_curvature,
    imaging,
    volume_n,
    what_age_was_the_patient_at_the_time_of_the_injury,
    did_the_patient_get_to_theatre_within_36_hours_of_presentation,
    what_age_was_the_patient_at_the_time_of_the_injury,
    was_the_patient_transferred_from_another_hospital_for_surgical_repair,
    was_a_urethral_injury_noted_at_surgery,
    laterality,
    distal,
    proximal,
    mid,
    crura,
    multi_level,
    mechanism_of_injury
  ) %>% drop_na()

abnormal_curvature_lr <- logistf(
  abnormal_curvature ~ did_the_patient_get_to_theatre_within_36_hours_of_presentation + imaging + volume_n + 
    what_age_was_the_patient_at_the_time_of_the_injury + was_a_urethral_injury_noted_at_surgery +
    was_the_patient_transferred_from_another_hospital_for_surgical_repair + distal + proximal + mid +
    crura + laterality + mechanism_of_injury,
  data = abnormal_curvature_for_log_reg_362
)



abnormal_curvature_lr %>%
  tbl_regression(
    exponentiate = TRUE,
    label = list(
      did_the_patient_get_to_theatre_within_36_hours_of_presentation = "Theatre within 36 hours of Presentation",
      was_a_urethral_injury_noted_at_surgery = "Urethral Injury",
      was_the_patient_transferred_from_another_hospital_for_surgical_repair = "Transfer for Surgical Repair",
      locality = "Site of Injury",
      distal = "Distal Corpus Location",
    proximal = "Proximal Corpus Location",
    mid = "Mid Corpus Location",
    crura = "Crural location",
      mechanism_of_injury = "Mechanism of Injury",
      what_age_was_the_patient_at_the_time_of_the_injury =
        "Age at Time of Injury",
    imaging = "Imaging",
    volume_n = "Centre Volume"
    )
  ) %>%
  bold_labels() %>%
  as_gt() %>%  
  tab_source_note(
    md(glue(
      "Number of patients in analysis: {nrow(abnormal_curvature_for_log_reg_362)}"
    ))
  )

## ! `broom::tidy()` failed to tidy the model.

## ✔ `tidy_parameters()` used instead.

## ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.

Characteristic	exp(Beta)	95% CI	p-value
Theatre within 36 hours of Presentation
No	—	—
Yes	0.49	0.17, 1.41	0.2
Imaging
No	—	—
Yes	0.36	0.11, 1.04	0.060
Centre Volume	1.03	0.99, 1.07	0.11
Age at Time of Injury	0.97	0.92, 1.01	0.11
Urethral Injury
No	—	—
Yes	1.07	0.35, 3.05	>0.9
Transfer for Surgical Repair
No	—	—
Yes	1.73	0.58, 5.25	0.3
Distal Corpus Location
No	—	—
Yes	0.39	0.00, 5.49	0.5
Proximal Corpus Location
No	—	—
Yes	0.39	0.00, 6.05	0.5
Mid Corpus Location
No	—	—
Yes	0.24	0.00, 3.35	0.3
Crural location
No	—	—
Yes	0.10	0.00, 8.89	0.3
laterality
Bilateral	—	—
Unilateral	0.57	0.17, 1.93	0.4
Mechanism of Injury
Sexual intercourse trauma	—	—
Other	0.29	0.00, 3.78	0.4
Non-sexual blunt force trauma	0.47	0.10, 1.92	0.3
Masturbation injury	0.12	0.00, 1.12	0.066
Abbreviation: CI = Confidence Interval
Number of patients in analysis: 128

8.4.3 ROC

abnormal_curvature_lr_predict <-
  predict(abnormal_curvature_lr, abnormal_curvature_for_log_reg_36, type = "response")

abnormal_curvature_outcome <- as.factor(ifelse(abnormal_curvature_for_log_reg_36$abnormal_curvature == "Yes",
                                  1,
                                  0))

abnormal_curvature_lr_roc <-
  pROC::roc(abnormal_curvature_outcome, abnormal_curvature_lr_predict)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

auc(abnormal_curvature_lr_roc)

## Area under the curve: 0.7712

ci.auc(abnormal_curvature_lr_roc)

## 95% CI: 0.6865-0.8558 (DeLong)

abnormal_curvature_lr_roc_data <- cbind(sensitivities = abnormal_curvature_lr_roc$sensitivities,
                                specificities = abnormal_curvature_lr_roc$specificities) %>% as_tibble()

ggplot(as_tibble(abnormal_curvature_lr_roc_data), aes(x = 1 - specificities, y = sensitivities)) +
  geom_line() +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed")

8.4.4 Calibration plot

abnormal_curvature_lr_predict2 <- as_tibble(abnormal_curvature_lr_predict) %>% cbind(abnormal_curvature_for_log_reg_36) %>% subset(select = c(
  value,
  abnormal_curvature,
  did_the_patient_get_to_theatre_within_36_hours_of_presentation
)) %>% mutate(
  pred = value,
  .keep = "unused"
)

abnormal_curvature_lr_predict2$abnormal_curvature <- ifelse(abnormal_curvature_lr_predict2$abnormal_curvature == "Yes", 1, 0)

abnormal_curvature_lr_predict2$did_the_patient_get_to_theatre_within_36_hours_of_presentation <- as.integer(abnormal_curvature_lr_predict2$did_the_patient_get_to_theatre_within_36_hours_of_presentation)

## Warning: NAs introduced by coercion

predtools::calibration_plot(
  data = abnormal_curvature_lr_predict2,
  obs = "abnormal_curvature",
  pred = "pred",
  y_lim = c(0, 1),
  x_lim = c(0, 1),
  title = "Calibration plot for abnormal_curvature Model"
)

## $calibration_plot

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

8.4.5 Brier Score

brier <- function(model, data, outcome) {
  probs <- predict(model, type = "response")
  obs <- data[[outcome]]
  mean((obs - probs)^2, na.rm = TRUE)
}

brier_score <- brier(abnormal_curvature_lr, abnormal_curvature_for_log_reg_362, "abnormal_curvature") 

cbind(
  "Score" = "Brier Score",
  "Brier Score" = round(brier_score, digits = 3)
) %>% as_tibble() %>% gt() %>% gt_theme_espn()

Score	Brier Score
Brier Score	0.169

8.4.6 Nagelkerke’s R squared

pseudoR2_logistf(abnormal_curvature_lr, abnormal_curvature_for_log_reg_362)

Score	Nagelkerke's R^2
Nagelkerke's R^2	0.238

9 Logistic Regression - Time from Presentation to Theatre with 48 hour cut-off

9.1 Pain

9.1.1 Post-hoc Power calculation

n1 <- nrow(pain_for_log_reg_48 %>% subset(did_the_patient_get_to_theatre_within_48_hours_of_presentation == "No"))
p1 <- nrow(pain_for_log_reg_48 %>% subset(did_the_patient_get_to_theatre_within_48_hours_of_presentation == "No" & pain == "Yes")) / n1

n2 <- nrow(pain_for_log_reg_48 %>% subset(did_the_patient_get_to_theatre_within_48_hours_of_presentation == "Yes"))
alpha <- 0.05

p2 <- nrow(pain_for_log_reg_48 %>% subset(did_the_patient_get_to_theatre_within_48_hours_of_presentation == "Yes" & pain == "Yes")) / n2

power <- calculate_power_table(p1, p2, n1, n2, alpha, title = "Power Calculation Results for Post-Operative Pain")
power

Power Calculation Results for Post-Operative Pain
Based on Two Proportion Test
Parameter	Value
Group 1 Proportion (p1)	0.03
Group 2 Proportion (p2)	0.13
Group 1 Size (n1)	32.00
Group 2 Size (n2)	213.00
Significance Level (alpha)	0.05
Calculated Power (%)	69.50

9.1.2 Logistic Regression (Firth Correction)

knife_to_skin1 <- knife_to_skin %>% 
  subset(select = c(
    submission_id,
    time_from_injury_to_th
  ))

pain_for_log_reg_48 <- pain_for_log_reg_48 %>% 
  left_join(knife_to_skin1,
            by = c("submission_id" = "submission_id")) %>%
  mutate(pain = pain %>% factor(levels = c("Yes",
                                           "No")))

pain_for_log_reg_482 <- pain_for_log_reg_48 %>%
  mutate(presentation_to_knife_to_skin = ifelse(
    presentation_to_knife_to_skin < 0 | presentation_to_knife_to_skin > 7920,
    NA,
    presentation_to_knife_to_skin / (60*6)
  ),
    pain = case_when(pain == "Yes" ~ 1,
                     pain == "No" ~ 0,
                     TRUE ~ NA_real_),
  .keep = "unused"
  ) %>%
  select(
    pain,
    imaging,
    volume_n,
    what_age_was_the_patient_at_the_time_of_the_injury,
    did_the_patient_get_to_theatre_within_48_hours_of_presentation,
    what_age_was_the_patient_at_the_time_of_the_injury,
    was_the_patient_transferred_from_another_hospital_for_surgical_repair,
    was_a_urethral_injury_noted_at_surgery,
    laterality,
    distal,
    proximal,
    mid,
    crura,
    multi_level,
    mechanism_of_injury
  ) %>% drop_na()

pain_lr <- logistf(
  pain ~ did_the_patient_get_to_theatre_within_48_hours_of_presentation + imaging + volume_n +
    what_age_was_the_patient_at_the_time_of_the_injury + was_a_urethral_injury_noted_at_surgery +
    was_the_patient_transferred_from_another_hospital_for_surgical_repair + distal + proximal + mid +
    crura + multi_level + laterality + mechanism_of_injury,
  data = pain_for_log_reg_482
)

pain_lr %>%
  tbl_regression(
    exponentiate = TRUE,
    label = list(
      did_the_patient_get_to_theatre_within_48_hours_of_presentation = "Theatre within 48 hours of Presentation",
      was_a_urethral_injury_noted_at_surgery = "Urethral Injury",
      was_the_patient_transferred_from_another_hospital_for_surgical_repair = "Transfer for Surgical Repair",
      laterality = "Side of Injury",
      distal = "Distal Corpus Location",
    proximal = "Proximal Corpus Location",
    mid = "Mid Corpus Location",
    crura = "Crural location",
    multi_level = "Multi-level Injury",
      mechanism_of_injury = "Mechanism of Injury",
      what_age_was_the_patient_at_the_time_of_the_injury =
        "Age at Time of Injury",
    imaging = "Imaging",
    volume_n = "Centre Volume"
    )
  ) %>%
  bold_labels() %>%
  as_gt() %>%  
  tab_source_note(
    md(glue(
      "Number of patients in analysis: {nrow(pain_for_log_reg_482)}"
    ))
  )

## ! `broom::tidy()` failed to tidy the model.

## ✔ `tidy_parameters()` used instead.

## ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.

Characteristic	exp(Beta)	95% CI	p-value
Theatre within 48 hours of Presentation
No	—	—
Yes	5.73	1.01, 86.9	0.048
Imaging
No	—	—
Yes	0.90	0.32, 2.60	0.8
Centre Volume	0.94	0.89, 0.98	0.005
Age at Time of Injury	1.04	0.99, 1.09	0.13
Urethral Injury
No	—	—
Yes	1.27	0.33, 4.32	0.7
Transfer for Surgical Repair
No	—	—
Yes	3.89	1.21, 13.5	0.023
Distal Corpus Location
No	—	—
Yes	0.03	0.00, 107	0.4
Proximal Corpus Location
No	—	—
Yes	0.04	0.00, 81.7	0.4
Mid Corpus Location
No	—	—
Yes	0.02	0.00, 43.0	0.3
Crural location
No	—	—
Yes	0.43	0.00, 263	0.8
Multi-level Injury
No	—	—
Yes	182	0.02, 3,307,433	0.3
Side of Injury
Bilateral	—	—
Unilateral	0.58	0.14, 2.66	0.5
Mechanism of Injury
Sexual intercourse trauma	—	—
Other	0.67	0.00, 8.98	0.8
Non-sexual blunt force trauma	3.80	0.77, 17.2	0.10
Masturbation injury	6.68	0.88, 60.1	0.066
Abbreviation: CI = Confidence Interval
Number of patients in analysis: 141

# Add back in to get 24hr cut-off: presentation_to_knife_to_skin,
# did_the_patient_get_to_theatre_within_48_hours_of_presentation
# time_from_injury_to_th

9.1.3 ROC

pain_lr_predict <-
  predict(pain_lr, pain_for_log_reg_48, type = "response")

pain_outcome <- as.factor(ifelse(pain_for_log_reg_48$pain == "Yes",
                                  1,
                                  0))

pain_lr_roc <-
  pROC::roc(pain_outcome, pain_lr_predict)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

auc(pain_lr_roc)

## Area under the curve: 0.8

ci.auc(pain_lr_roc)

## 95% CI: 0.7044-0.8957 (DeLong)

pain_lr_roc_data <- cbind(sensitivities = pain_lr_roc$sensitivities,
                                specificities = pain_lr_roc$specificities) %>% as_tibble()

ggplot(as_tibble(pain_lr_roc_data), aes(x = 1 - specificities, y = sensitivities)) +
  geom_line() +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed")

9.1.4 Calibration plot

pain_lr_predict2 <- as_tibble(pain_lr_predict) %>% cbind(pain_for_log_reg_48) %>% subset(select = c(
  value,
  pain,
  did_the_patient_get_to_theatre_within_48_hours_of_presentation
)) %>% mutate(
  pred = value,
  .keep = "unused"
)

pain_lr_predict2$pain <- ifelse(pain_lr_predict2$pain == "Yes", 1, 0)

pain_lr_predict2$did_the_patient_get_to_theatre_within_48_hours_of_presentation <- as.integer(pain_lr_predict2$did_the_patient_get_to_theatre_within_48_hours_of_presentation)

## Warning: NAs introduced by coercion

predtools::calibration_plot(
  data = pain_lr_predict2,
  obs = "pain",
  pred = "pred",
  y_lim = c(0, 1),
  x_lim = c(0, 1),
  title = "Calibration plot for Pain Model"
)

## $calibration_plot

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

9.1.5 Brier Score

brier <- function(model, data, outcome) {
  probs <- predict(model, type = "response")
  obs <- data[[outcome]]
  mean((obs - probs)^2, na.rm = TRUE)
}

brier_score <- brier(pain_lr, pain_for_log_reg_482, "pain") 

cbind(
  "Score" = "Brier Score",
  "Brier Score" = round(brier_score, digits = 3)
) %>% as_tibble() %>% gt() %>% gt_theme_espn()

Score	Brier Score
Brier Score	0.112

9.1.6 Nagelkerke’s R squared

pseudoR2_logistf <- function(model, data) {
  # Extract log-likelihoods
  ll_full <- model$loglik["full"]
  ll_null <- model$loglik["null"]
  
  # Sample size
  n <- nrow(data)
  
  # Cox & Snell R²
  R2_CS <- 1 - exp((2 / n) * (ll_null - ll_full))
  
  # Nagelkerke R²
  R2_N <- R2_CS / (1 - exp((2 / n) * ll_null))

  R2_N <- R2_N %>% as_tibble()
  
  x <- cbind(
  "Score" = "Nagelkerke's R^2",
  "Nagelkerke's R^2" = round(R2_N$value, digits = 3)
) %>% as_tibble() %>% gt() %>% gt_theme_espn()
  
  # Return as a named list
  return(
    x
  )
}

pseudoR2_logistf(pain_lr, pain_for_log_reg_482)

Score	Nagelkerke's R^2
Nagelkerke's R^2	0.288

9.2 ED

9.2.1 Post-hoc Power calculation

n1 <- nrow(ed_for_log_reg_48 %>% subset(did_the_patient_get_to_theatre_within_48_hours_of_presentation == "No"))
p1 <- nrow(ed_for_log_reg_48 %>% subset(did_the_patient_get_to_theatre_within_48_hours_of_presentation == "No" & ed == "Yes")) / n1

n2 <- nrow(ed_for_log_reg_48 %>% subset(did_the_patient_get_to_theatre_within_48_hours_of_presentation == "Yes"))
alpha <- 0.05

p2 <- nrow(ed_for_log_reg_48 %>% subset(did_the_patient_get_to_theatre_within_48_hours_of_presentation == "Yes" & ed == "Yes")) / n2

power <- calculate_power_table(p1, p2, n1, n2, alpha, title = "Power Calculation Results for Post-Operative Erectile Dysfunction")
power

Power Calculation Results for Post-Operative Erectile Dysfunction
Based on Two Proportion Test
Parameter	Value
Group 1 Proportion (p1)	0.28
Group 2 Proportion (p2)	0.14
Group 1 Size (n1)	32.00
Group 2 Size (n2)	213.00
Significance Level (alpha)	0.05
Calculated Power (%)	0.00

9.2.2 Logistic Regression (Firth Correction)

ed_for_log_reg_48 <- ed_for_log_reg_48 %>%
  left_join(knife_to_skin1, by = c("submission_id" = "submission_id")) %>%
  mutate(ed = ed %>% factor(levels = c("No", "Yes")))

ed_for_log_reg_482 <- ed_for_log_reg_48 %>%
  drop_na(
    ed,
    imaging,
    volume_n,
    did_the_patient_get_to_theatre_within_48_hours_of_presentation,
    what_age_was_the_patient_at_the_time_of_the_injury,
    presentation_to_knife_to_skin,
    what_age_was_the_patient_at_the_time_of_the_injury,
    was_the_patient_transferred_from_another_hospital_for_surgical_repair,
    was_a_urethral_injury_noted_at_surgery,
    laterality,
    distal,
    proximal,
    mid,
    crura,
    multi_level,
    mechanism_of_injury
  ) %>%
  mutate(
    ed = case_when(ed == "Yes" ~ 1,
                   ed == "No" ~ 0,
                   TRUE ~ NA_real_),
    presentation_to_knife_to_skin = ifelse(
    presentation_to_knife_to_skin < 0 | presentation_to_knife_to_skin > 5000,
    NA,
    presentation_to_knife_to_skin / (60*6)
  ),
    .keep = "unused"
  )

ed_lr <- logistf(
  ed ~ did_the_patient_get_to_theatre_within_48_hours_of_presentation + imaging + volume_n +
    what_age_was_the_patient_at_the_time_of_the_injury + was_a_urethral_injury_noted_at_surgery +
    was_the_patient_transferred_from_another_hospital_for_surgical_repair + distal + proximal + mid +
    crura + multi_level + laterality + mechanism_of_injury,
  data = ed_for_log_reg_482
)

ed_lr %>%
  tbl_regression(
    exponentiate = TRUE,
    label = list(
      did_the_patient_get_to_theatre_within_48_hours_of_presentation = "Theatre within 48 hours of Presentation",
      was_a_urethral_injury_noted_at_surgery = "Urethral Injury",
      was_the_patient_transferred_from_another_hospital_for_surgical_repair = "Transfer for Surgical Repair",
      distal = "Distal Corpus Location",
    proximal = "Proximal Corpus Location",
    mid = "Mid Corpus Location",
    crura = "Crural location",
    multi_level = "Multi-level Injury",
      laterality = "Side of Injury",
      mechanism_of_injury = "Mechanism of Injury",
      what_age_was_the_patient_at_the_time_of_the_injury =
        "Age at Time of Injury",
    imaging = "Imaging",
    volume_n = "Centre Volume"
    )
  ) %>%
  bold_labels() %>%
  as_gt() %>%  
  tab_source_note(
    md(glue(
      "Number of patients in analysis: {nrow(ed_for_log_reg_482)}"
    ))
  )

## ! `broom::tidy()` failed to tidy the model.

## ✔ `tidy_parameters()` used instead.

## ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.

Characteristic	exp(Beta)	95% CI	p-value
Theatre within 48 hours of Presentation
No	—	—
Yes	0.30	0.09, 0.98	0.047
Imaging
No	—	—
Yes	1.22	0.45, 3.49	0.7
Centre Volume	0.98	0.95, 1.02	0.3
Age at Time of Injury	1.02	0.98, 1.07	0.3
Urethral Injury
No	—	—
Yes	2.97	0.98, 9.44	0.054
Transfer for Surgical Repair
No	—	—
Yes	1.82	0.58, 6.10	0.3
Distal Corpus Location
No	—	—
Yes	19.6	0.02, 58,247	0.4
Proximal Corpus Location
No	—	—
Yes	15.9	0.03, 36,660	0.4
Mid Corpus Location
No	—	—
Yes	6.79	0.01, 15,843	0.6
Crural location
No	—	—
Yes	7.39	0.03, 6,895	0.5
Multi-level Injury
No	—	—
Yes	0.14	0.00, 380	0.6
Side of Injury
Bilateral	—	—
Unilateral	6.24	1.35, 45.5	0.017
Mechanism of Injury
Sexual intercourse trauma	—	—
Other	0.30	0.00, 4.68	0.4
Non-sexual blunt force trauma	0.29	0.03, 1.38	0.13
Masturbation injury	0.11	0.00, 1.20	0.076
Abbreviation: CI = Confidence Interval
Number of patients in analysis: 136

# did_the_patient_get_to_theatre_within_48_hours_of_presentation,

9.2.3 ROC

ed_lr_predict <-
  predict(ed_lr, ed_for_log_reg_48, type = "response")

ed_outcome <- as.factor(ifelse(ed_for_log_reg_48$ed == "Yes",
                                  1,
                                  0))

ed_lr_roc <-
  pROC::roc(ed_outcome, ed_lr_predict)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

auc(ed_lr_roc)

## Area under the curve: 0.8157

ci.auc(ed_lr_roc)

## 95% CI: 0.7381-0.8933 (DeLong)

ed_lr_roc_data <- cbind(sensitivities = ed_lr_roc$sensitivities,
                                specificities = ed_lr_roc$specificities) %>% as_tibble()

ggplot(as_tibble(ed_lr_roc_data), aes(x = 1 - specificities, y = sensitivities)) +
  geom_line() +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed")

9.2.4 Calibration plot

ed_lr_predict2 <- as_tibble(ed_lr_predict) %>% cbind(ed_for_log_reg_48) %>% subset(select = c(
  value,
  ed,
  did_the_patient_get_to_theatre_within_48_hours_of_presentation
)) %>% mutate(
  pred = value,
  .keep = "unused"
)

ed_lr_predict2$ed <- ifelse(ed_lr_predict2$ed == "Yes", 1, 0)

ed_lr_predict2$did_the_patient_get_to_theatre_within_48_hours_of_presentation <- as.integer(ed_lr_predict2$did_the_patient_get_to_theatre_within_48_hours_of_presentation)

## Warning: NAs introduced by coercion

predtools::calibration_plot(
  data = ed_lr_predict2,
  obs = "ed",
  pred = "pred",
  y_lim = c(0, 1),
  x_lim = c(0, 1),
  title = "Calibration plot for ed Model"
)

## $calibration_plot

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

9.2.5 Brier Score

brier_score <- brier(ed_lr, ed_for_log_reg_482, "ed")

cbind(
  "Score" = "Brier Score",
  "Brier Score" = round(brier_score, digits = 3)
) %>% as_tibble() %>% gt() %>% gt_theme_espn()

Score	Brier Score
Brier Score	0.154

9.2.6 Nagelkerke’s R squared

pseudoR2_logistf(ed_lr, ed_for_log_reg_482)

Score	Nagelkerke's R^2
Nagelkerke's R^2	0.268

9.3 Waisting / Shortening

9.3.1 Post-hoc Power calculation

n1 <- nrow(waisting_shortening_for_log_reg_48 %>% subset(did_the_patient_get_to_theatre_within_48_hours_of_presentation == "No"))
p1 <- nrow(waisting_shortening_for_log_reg_48 %>% subset(did_the_patient_get_to_theatre_within_48_hours_of_presentation == "No" & waisting_shortening == "Yes")) / n1

n2 <- nrow(waisting_shortening_for_log_reg_48 %>% subset(did_the_patient_get_to_theatre_within_48_hours_of_presentation == "Yes"))
alpha <- 0.05

p2 <- nrow(waisting_shortening_for_log_reg_48 %>% subset(did_the_patient_get_to_theatre_within_48_hours_of_presentation == "Yes" & waisting_shortening == "Yes")) / n2

power <- calculate_power_table(p1, p2, n1, n2, alpha, title = "Power Calculation Results for Post-Operative waisting_shortening")
power

Power Calculation Results for Post-Operative waisting_shortening
Based on Two Proportion Test
Parameter	Value
Group 1 Proportion (p1)	0.09
Group 2 Proportion (p2)	0.07
Group 1 Size (n1)	32.00
Group 2 Size (n2)	213.00
Significance Level (alpha)	0.05
Calculated Power (%)	0.70

9.3.2 Logistic Regression (Firth Correction)

waisting_shortening_for_log_reg_48 <- waisting_shortening_for_log_reg_48 %>% 
  left_join(knife_to_skin1,
            by = c("submission_id" = "submission_id")) %>%
  mutate(waisting_shortening = waisting_shortening %>% factor(levels = c("Yes",
                                           "No")))

waisting_shortening_for_log_reg_482 <- waisting_shortening_for_log_reg_48 %>%
  mutate(presentation_to_knife_to_skin = ifelse(
    presentation_to_knife_to_skin < 0 | presentation_to_knife_to_skin > 5000,
    NA,
    presentation_to_knife_to_skin / (60*6)
  ),
    waisting_shortening = case_when(waisting_shortening == "Yes" ~ 1,
                     waisting_shortening == "No" ~ 0,
                     TRUE ~ NA_real_),
  .keep = "unused"
  ) %>%
  select(
    waisting_shortening,
    imaging,
    volume_n,
    what_age_was_the_patient_at_the_time_of_the_injury,
    did_the_patient_get_to_theatre_within_48_hours_of_presentation,
    what_age_was_the_patient_at_the_time_of_the_injury,
    was_the_patient_transferred_from_another_hospital_for_surgical_repair,
    was_a_urethral_injury_noted_at_surgery,
    laterality,
    locality,
    mechanism_of_injury
  ) %>% drop_na()

waisting_shortening_lr <- logistf(
  waisting_shortening ~ did_the_patient_get_to_theatre_within_48_hours_of_presentation + imaging + volume_n +
    what_age_was_the_patient_at_the_time_of_the_injury + was_a_urethral_injury_noted_at_surgery +
    was_the_patient_transferred_from_another_hospital_for_surgical_repair + locality + laterality + mechanism_of_injury,
  data = waisting_shortening_for_log_reg_482
)

waisting_shortening_lr %>%
  tbl_regression(
    exponentiate = TRUE,
    label = list(
      did_the_patient_get_to_theatre_within_48_hours_of_presentation = "Theatre within 48 hours of Presentation",
      was_a_urethral_injury_noted_at_surgery = "Urethral Injury",
      was_the_patient_transferred_from_another_hospital_for_surgical_repair = "Transfer for Surgical Repair",
      locality = "Site of Injury",
      laterality = "Side of Injury",
      mechanism_of_injury = "Mechanism of Injury",
      what_age_was_the_patient_at_the_time_of_the_injury =
        "Age at Time of Injury",
      imaging = "Imaging",
      volume_n = "Centre Volume"
    )
  ) %>%
  bold_labels() %>%
  as_gt() %>%  
  tab_source_note(
    md(glue(
      "Number of patients in analysis: {nrow(waisting_shortening_for_log_reg_482)}"
    ))
  )

## ! `broom::tidy()` failed to tidy the model.

## ✔ `tidy_parameters()` used instead.

## ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.

Characteristic	exp(Beta)	95% CI	p-value
Theatre within 48 hours of Presentation
No	—	—
Yes	0.52	0.12, 2.47	0.4
Imaging
No	—	—
Yes	0.49	0.12, 1.86	0.3
Centre Volume	1.02	0.97, 1.07	0.5
Age at Time of Injury	1.00	0.95, 1.06	0.9
Urethral Injury
No	—	—
Yes	0.95	0.18, 3.85	>0.9
Transfer for Surgical Repair
No	—	—
Yes	0.62	0.12, 2.94	0.6
Site of Injury
Crura	—	—
Distal	0.39	0.01, 84.5	0.7
Mid	0.76	0.02, 145	0.9
Multi-level	0.47	0.00, 130	0.8
Proximal	0.69	0.02, 130	0.8
Side of Injury
Bilateral	—	—
Unilateral	0.31	0.07, 1.32	0.11
Mechanism of Injury
Sexual intercourse trauma	—	—
Other	0.61	0.00, 8.60	0.7
Non-sexual blunt force trauma	0.23	0.00, 2.18	0.2
Masturbation injury	0.21	0.00, 2.93	0.3
Abbreviation: CI = Confidence Interval
Number of patients in analysis: 119

9.3.3 ROC

waisting_shortening_lr_predict <-
  predict(waisting_shortening_lr, waisting_shortening_for_log_reg_48, type = "response")

waisting_shortening_outcome <- as.factor(ifelse(waisting_shortening_for_log_reg_48$waisting_shortening == "Yes",
                                  1,
                                  0))

waisting_shortening_lr_roc <-
  pROC::roc(waisting_shortening_outcome, waisting_shortening_lr_predict)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

auc(waisting_shortening_lr_roc)

## Area under the curve: 0.7699

ci.auc(waisting_shortening_lr_roc)

## 95% CI: 0.6546-0.8852 (DeLong)

waisting_shortening_lr_roc_data <- cbind(sensitivities = waisting_shortening_lr_roc$sensitivities,
                                specificities = waisting_shortening_lr_roc$specificities) %>% as_tibble()

ggplot(as_tibble(waisting_shortening_lr_roc_data), aes(x = 1 - specificities, y = sensitivities)) +
  geom_line() +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed")

9.3.4 Calibration plot

waisting_shortening_lr_predict2 <- as_tibble(waisting_shortening_lr_predict) %>% cbind(waisting_shortening_for_log_reg_48) %>% subset(select = c(
  value,
  waisting_shortening,
  did_the_patient_get_to_theatre_within_48_hours_of_presentation
)) %>% mutate(
  pred = value,
  .keep = "unused"
)

waisting_shortening_lr_predict2$waisting_shortening <- ifelse(waisting_shortening_lr_predict2$waisting_shortening == "Yes", 1, 0)

waisting_shortening_lr_predict2$did_the_patient_get_to_theatre_within_48_hours_of_presentation <- as.integer(waisting_shortening_lr_predict2$did_the_patient_get_to_theatre_within_48_hours_of_presentation)

## Warning: NAs introduced by coercion

predtools::calibration_plot(
  data = waisting_shortening_lr_predict2,
  obs = "waisting_shortening",
  pred = "pred",
  y_lim = c(0, 1),
  x_lim = c(0, 1),
  title = "Calibration plot for waisting_shortening Model"
)

## $calibration_plot

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

9.3.5 Brier Score

brier <- function(model, data, outcome) {
  probs <- predict(model, type = "response")
  obs <- data[[outcome]]
  mean((obs - probs)^2, na.rm = TRUE)
}

brier_score <- brier(waisting_shortening_lr, waisting_shortening_for_log_reg_482, "waisting_shortening") 

cbind(
  "Score" = "Brier Score",
  "Brier Score" = round(brier_score, digits = 3)
) %>% as_tibble() %>% gt() %>% gt_theme_espn()

Score	Brier Score
Brier Score	0.102

9.3.6 Nagelkerke’s R squared

pseudoR2_logistf(waisting_shortening_lr, waisting_shortening_for_log_reg_482)

Score	Nagelkerke's R^2
Nagelkerke's R^2	0.14

9.4 Abnormal Curvature

9.4.1 Post-hoc Power calculation

n1 <- nrow(abnormal_curvature_for_log_reg_48 %>% subset(did_the_patient_get_to_theatre_within_48_hours_of_presentation == "No"))
p1 <- nrow(abnormal_curvature_for_log_reg_48 %>% subset(did_the_patient_get_to_theatre_within_48_hours_of_presentation == "No" & abnormal_curvature == "Yes")) / n1

n2 <- nrow(abnormal_curvature_for_log_reg_48 %>% subset(did_the_patient_get_to_theatre_within_48_hours_of_presentation == "Yes"))
alpha <- 0.05

p2 <- nrow(abnormal_curvature_for_log_reg_48 %>% subset(did_the_patient_get_to_theatre_within_48_hours_of_presentation == "Yes" & abnormal_curvature == "Yes")) / n2

power <- calculate_power_table(p1, p2, n1, n2, alpha, title = "Power Calculation Results for Post-Operative abnormal_curvature")
power

Power Calculation Results for Post-Operative abnormal_curvature
Based on Two Proportion Test
Parameter	Value
Group 1 Proportion (p1)	0.22
Group 2 Proportion (p2)	0.15
Group 1 Size (n1)	32.00
Group 2 Size (n2)	213.00
Significance Level (alpha)	0.05
Calculated Power (%)	0.20

9.4.2 Logistic Regression (Firth Correction)

Multi-level removed as model failed to separate

abnormal_curvature_for_log_reg_48 <- abnormal_curvature_for_log_reg_48 %>% 
  left_join(knife_to_skin1,
            by = c("submission_id" = "submission_id")) %>%
  mutate(abnormal_curvature = abnormal_curvature %>% factor(levels = c("Yes",
                                           "No")))

abnormal_curvature_for_log_reg_482 <- abnormal_curvature_for_log_reg_48 %>%
  mutate(presentation_to_knife_to_skin = ifelse(
    presentation_to_knife_to_skin < 0 | presentation_to_knife_to_skin > 5000,
    NA,
    presentation_to_knife_to_skin / (60*6)
  ),
    abnormal_curvature = case_when(abnormal_curvature == "Yes" ~ 1,
                     abnormal_curvature == "No" ~ 0,
                     TRUE ~ NA_real_),
  .keep = "unused"
  ) %>%
  select(
    abnormal_curvature,
    imaging,
    volume_n,
    what_age_was_the_patient_at_the_time_of_the_injury,
    did_the_patient_get_to_theatre_within_48_hours_of_presentation,
    what_age_was_the_patient_at_the_time_of_the_injury,
    was_the_patient_transferred_from_another_hospital_for_surgical_repair,
    was_a_urethral_injury_noted_at_surgery,
    laterality,
    distal,
    proximal,
    mid,
    crura,
    multi_level,
    mechanism_of_injury
  ) %>% drop_na()

abnormal_curvature_lr <- logistf(
  abnormal_curvature ~ did_the_patient_get_to_theatre_within_48_hours_of_presentation + imaging + volume_n + 
    what_age_was_the_patient_at_the_time_of_the_injury + was_a_urethral_injury_noted_at_surgery +
    was_the_patient_transferred_from_another_hospital_for_surgical_repair + distal + proximal + mid +
    crura + laterality + mechanism_of_injury,
  data = abnormal_curvature_for_log_reg_482
)



abnormal_curvature_lr %>%
  tbl_regression(
    exponentiate = TRUE,
    label = list(
      did_the_patient_get_to_theatre_within_48_hours_of_presentation = "Theatre within 48 hours of Presentation",
      was_a_urethral_injury_noted_at_surgery = "Urethral Injury",
      was_the_patient_transferred_from_another_hospital_for_surgical_repair = "Transfer for Surgical Repair",
      locality = "Site of Injury",
      distal = "Distal Corpus Location",
    proximal = "Proximal Corpus Location",
    mid = "Mid Corpus Location",
    crura = "Crural location",
      mechanism_of_injury = "Mechanism of Injury",
      what_age_was_the_patient_at_the_time_of_the_injury =
        "Age at Time of Injury",
    imaging = "Imaging",
    volume_n = "Centre Volume"
    )
  ) %>%
  bold_labels() %>%
  as_gt() %>%  
  tab_source_note(
    md(glue(
      "Number of patients in analysis: {nrow(abnormal_curvature_for_log_reg_482)}"
    ))
  )

## ! `broom::tidy()` failed to tidy the model.

## ✔ `tidy_parameters()` used instead.

## ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.

Characteristic	exp(Beta)	95% CI	p-value
Theatre within 48 hours of Presentation
No	—	—
Yes	0.39	0.12, 1.23	0.11
Imaging
No	—	—
Yes	0.38	0.13, 1.06	0.066
Centre Volume	1.03	0.99, 1.07	0.11
Age at Time of Injury	0.97	0.93, 1.01	0.14
Urethral Injury
No	—	—
Yes	1.09	0.36, 3.12	0.9
Transfer for Surgical Repair
No	—	—
Yes	1.75	0.58, 5.37	0.3
Distal Corpus Location
No	—	—
Yes	0.38	0.00, 5.28	0.5
Proximal Corpus Location
No	—	—
Yes	0.42	0.00, 6.69	0.6
Mid Corpus Location
No	—	—
Yes	0.27	0.00, 3.79	0.4
Crural location
No	—	—
Yes	0.11	0.00, 10.2	0.4
laterality
Bilateral	—	—
Unilateral	0.56	0.17, 1.87	0.3
Mechanism of Injury
Sexual intercourse trauma	—	—
Other	0.32	0.00, 3.96	0.4
Non-sexual blunt force trauma	0.47	0.10, 1.90	0.3
Masturbation injury	0.10	0.00, 1.06	0.057
Abbreviation: CI = Confidence Interval
Number of patients in analysis: 128

9.4.3 ROC

abnormal_curvature_lr_predict <-
  predict(abnormal_curvature_lr, abnormal_curvature_for_log_reg_48, type = "response")

abnormal_curvature_outcome <- as.factor(ifelse(abnormal_curvature_for_log_reg_48$abnormal_curvature == "Yes",
                                  1,
                                  0))

abnormal_curvature_lr_roc <-
  pROC::roc(abnormal_curvature_outcome, abnormal_curvature_lr_predict)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

auc(abnormal_curvature_lr_roc)

## Area under the curve: 0.778

ci.auc(abnormal_curvature_lr_roc)

## 95% CI: 0.696-0.8599 (DeLong)

abnormal_curvature_lr_roc_data <- cbind(sensitivities = abnormal_curvature_lr_roc$sensitivities,
                                specificities = abnormal_curvature_lr_roc$specificities) %>% as_tibble()

ggplot(as_tibble(abnormal_curvature_lr_roc_data), aes(x = 1 - specificities, y = sensitivities)) +
  geom_line() +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed")

9.4.4 Calibration plot

abnormal_curvature_lr_predict2 <- as_tibble(abnormal_curvature_lr_predict) %>% cbind(abnormal_curvature_for_log_reg_48) %>% subset(select = c(
  value,
  abnormal_curvature,
  did_the_patient_get_to_theatre_within_48_hours_of_presentation
)) %>% mutate(
  pred = value,
  .keep = "unused"
)

abnormal_curvature_lr_predict2$abnormal_curvature <- ifelse(abnormal_curvature_lr_predict2$abnormal_curvature == "Yes", 1, 0)

abnormal_curvature_lr_predict2$did_the_patient_get_to_theatre_within_48_hours_of_presentation <- as.integer(abnormal_curvature_lr_predict2$did_the_patient_get_to_theatre_within_48_hours_of_presentation)

## Warning: NAs introduced by coercion

predtools::calibration_plot(
  data = abnormal_curvature_lr_predict2,
  obs = "abnormal_curvature",
  pred = "pred",
  y_lim = c(0, 1),
  x_lim = c(0, 1),
  title = "Calibration plot for abnormal_curvature Model"
)

## $calibration_plot

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

9.4.5 Brier Score

brier <- function(model, data, outcome) {
  probs <- predict(model, type = "response")
  obs <- data[[outcome]]
  mean((obs - probs)^2, na.rm = TRUE)
}

brier_score <- brier(abnormal_curvature_lr, abnormal_curvature_for_log_reg_482, "abnormal_curvature") 

cbind(
  "Score" = "Brier Score",
  "Brier Score" = round(brier_score, digits = 3)
) %>% as_tibble() %>% gt() %>% gt_theme_espn()

Score	Brier Score
Brier Score	0.168

9.4.6 Nagelkerke’s R squared

pseudoR2_logistf(abnormal_curvature_lr, abnormal_curvature_for_log_reg_482)

Score	Nagelkerke's R^2
Nagelkerke's R^2	0.247