SNAP Audit
R Geraghty
2025-09-12
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))## 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 Pain
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_
))
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.2 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.3 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.4 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.5 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.6 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.7 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)) ~ "Left",
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)) ~ "Right",
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 798 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 797 remaining warnings.2.8 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 Age
3.2.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.2.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.3 Type of Suture
Could make this into composite of suture type or suture diameter but
unusable as factor in current format
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 2-0 Suture | 93 | 38 |
| PDS 3-0 Suture | 46 | 19 |
| Other | 30 | 12 |
| Vicryl 2-0 Suture | 10 | 4 |
| Vicryl 3-0 Suture | 5 | 2 |
| PDS 0-Suture - PDS 2-0 Suture | 2 | 1 |
| PDS 0-Suture | 37 | 15 |
| PDS 1-0 Suture | 4 | 2 |
| Vicryl 0 -Suture | 2 | 1 |
| PDS 2-0 Suture - Vicryl 3-0 Suture | 1 | 0 |
| Vicryl 4-0 Suture | 3 | 1 |
| PDS 4-0 Suture | 10 | 4 |
| PDS 2-0 Suture - Vicryl 2-0 Suture | 1 | 0 |
| Nylon 4-0 Suture | 1 | 0 |
| Nylon 3-0 Suture | 1 | 0 |
3.4 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.5 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.6 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.7 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.8 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 |
|---|---|---|
| Yes | 191 | 85 |
| No | 35 | 15 |
3.9 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 |
|---|---|---|
| No | 25 | 12 |
| Yes | 179 | 88 |
3.10 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.11 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)| 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.12 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 |
|---|---|---|
| Subcoronal with degloving of penis | 81 | 32 |
| Penoscrotal | 96 | 38 |
| longitudinal | 59 | 24 |
| Subcoronal without degloving of penis | 7 | 3 |
| Other | 8 | 3 |
3.13 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.14 Time Variables
3.14.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 61 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.14.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 47 rows containing non-finite outside the scale range
## (`stat_bin()`).
3.14.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.15 Imaging
3.15.1 Correlation between Imaging Performed and Reporting Radiologist Subspecialism
3.15.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.15.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.15.2 Correlation between Imaging Performed and Transfer for Imaging
3.15.2.1 Any 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) %>%
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.15.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.15.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.15.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.15.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.16 Volume
3.16.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.16.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$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]", " "
))) %>%
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_
)
)
predictors$volume <- factor(predictors$volume,
levels = c("High",
"Medium",
"Low"))
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.17 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", "Left", "Right", "Total")) %>%
kable_styling(full_width = FALSE,
bootstrap_options = c("striped", "hover", "condensed")) %>%
row_spec(0, bold = TRUE, background = "#D3D3D3")| Location | Bilateral | Left | Right | Total |
|---|---|---|---|---|
| Crura | 0 (0%) | 4 (44.4%) | 5 (55.6%) | 9 (100%) |
| Distal | 1 (4.3%) | 9 (39.1%) | 13 (56.5%) | 23 (100%) |
| Mid | 7 (9.5%) | 27 (36.5%) | 40 (54.1%) | 74 (100%) |
| Multi-level | 3 (42.9%) | 2 (28.6%) | 2 (28.6%) | 7 (100%) |
| Proximal | 14 (12.3%) | 40 (35.1%) | 60 (52.6%) | 114 (100%) |
| Total | 25 (11%) | 82 (36.1%) | 120 (52.9%) | 227 (100%) |
3.18 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 | hospital_clean | volume | proximal | mid | distal | crura | multi_level | imaging |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 6230380745395251275 | 49 | Sexual intercourse trauma | NA | No | Yes | No | No | Yes | Scrotum | No | No clinical suspicion of urethral injury | No | No | NA | Yes | Yes | No | NA | Yes | Yes | Consultant | NA | No | NA | Subcoronal with degloving of penis | NA | Not performed | No | NA | Yes | 840 | NA | NA | West Midlands Alexandra Hospital, Redditch | Low | No | No | No | No | No | No |
| 6225338547719659855 | 56 | Sexual intercourse trauma | NA | No | No | No | Yes | Yes | Lower abdomen | No | No clinical suspicion of urethral injury | No | No | NA | NA | NA | No | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | Yes | NA | NA | NA | East Midlands Nottingham City Hospital | Medium | No | No | No | No | No | No |
| 6225308917713028889 | 34 | Sexual intercourse trauma | NA | No | NA | Yes | NA | Yes | NA | No | No clinical suspicion of urethral injury | No | No | NA | Yes | Yes | No | Yes | Yes | Yes | Consultant | Yes | Yes | Proximal shaft right corpora | Penoscrotal | No | Flexible cystoscopy | No | NA | Yes | 180 | Right | Proximal | East Midlands Nottingham City Hospital | Medium | Yes | No | No | No | No | No |
| 6225300887711789701 | 66 | Sexual intercourse trauma | Patient on top (Missionary/face to face) | No | No | Yes | NA | Yes | Scrotum - Lower abdomen | No | No clinical suspicion of urethral injury | No | NA | NA | No | Yes | Yes | Yes | Yes | No | Consultant | Yes | Yes | Proximal shaft left corpora | longitudinal | No | Not performed | No | NA | Yes | 635 | Left | Proximal | East Midlands Nottingham City Hospital | Medium | Yes | No | No | No | No | No |
| 6225291667712045628 | 56 | NA | NA | No | NA | NA | NA | NA | NA | No | NA | No | No | NA | No | Yes | Yes | Yes | Yes | Yes | Consultant | Yes | No | NA | Subcoronal with degloving of penis | No | Flexible cystoscopy | No | NA | Yes | 55 | NA | NA | East Midlands Nottingham City Hospital | Medium | No | No | No | No | No | No |
| 6225154497717063495 | 71 | Sexual intercourse trauma | NA | No | Yes | Yes | No | Yes | NA | No | No clinical suspicion of urethral injury | MRI | No | Yes | No | No | No | Yes | Yes | Yes | Consultant - Specialist Registrar | Yes | Yes | Distal shaft right corpora | Subcoronal without degloving of penis | Yes | Flexible cystoscopy | No | NA | Yes | 7310 | Right | Distal | East Midlands Nottingham City Hospital | Medium | No | No | Yes | No | No | Yes |
3.19 Plot missing
plot_missing(predictors_pres_to_th)
3.20 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
))
plot_missing(pain_for_log_reg)
glue("Number with complete data: {nrow(pain_for_log_reg %>% drop_na())}")## Number with complete data: 564.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
))
plot_missing(ed_for_log_reg)
glue("Number with complete data: {nrow(ed_for_log_reg %>% drop_na())}")## Number with complete data: 534.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
))
plot_missing(waisting_shortening_for_log_reg)
glue("Number with complete data: {nrow(waisting_shortening_for_log_reg %>% drop_na())}")## Number with complete data: 524.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
))
plot_missing(abnormal_curvature_for_log_reg)
glue("Number with complete data: {nrow(abnormal_curvature_for_log_reg %>% drop_na())}")## Number with complete data: 514.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 25 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 25 rows containing non-finite outside the scale range
## (`stat_density()`).## Warning: Removed 6 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, fill = ed)) +
geom_histogram(bins = 50) +
labs(
title = "Time from Presentation to Theatre",
x = "Minutes",
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, fill = ed)) +
geom_histogram(bins = 50) +
labs(
title = "Time from Injury to Theatre",
x = "Minutes",
y = "Count"
) +
theme_minimal() + xlim(0,100)
cowplot::plot_grid(p1, p2)## Warning: Removed 143 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 127 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()`).
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, fill = ed)) +
geom_density() +
labs(
title = "Time from Presentation to Theatre",
x = "Minutes",
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, fill = ed)) +
geom_density() +
labs(
title = "Time from Injury to Theatre",
x = "Minutes",
y = "Count"
) +
theme_minimal() + xlim(0,100)
cowplot::plot_grid(p1, p2)## Warning: Removed 143 rows containing non-finite outside the scale range
## (`stat_density()`).## Warning: Removed 127 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 18 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 18 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 24 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 24 rows containing non-finite outside the scale range
## (`stat_density()`).## Warning: Removed 5 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,
volume,
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 +
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 = "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 |
|---|---|---|---|
| Theatre within 24 hours of Presentation | |||
| Yes | — | — | |
| No | 0.81 | 0.24, 2.56 | 0.7 |
| Imaging | |||
| No | — | — | |
| Yes | 1.08 | 0.32, 3.86 | >0.9 |
| Centre Volume | |||
| High | — | — | |
| Medium | 8.94 | 1.23, 126 | 0.030 |
| Low | 11.9 | 2.09, 149 | 0.004 |
| Age at Time of Injury | 1.04 | 0.99, 1.10 | 0.14 |
| Urethral Injury | |||
| No | — | — | |
| Yes | 1.00 | 0.24, 3.55 | >0.9 |
| Transfer for Surgical Repair | |||
| No | — | — | |
| Yes | 6.06 | 1.64, 34.5 | 0.006 |
| Distal Corpus Location | |||
| No | — | — | |
| Yes | 0.00 | 0.00, 16.1 | 0.2 |
| Proximal Corpus Location | |||
| No | — | — | |
| Yes | 0.00 | 0.00, 11.7 | 0.2 |
| Mid Corpus Location | |||
| No | — | — | |
| Yes | 0.00 | 0.00, 6.18 | 0.14 |
| Crural location | |||
| No | — | — | |
| Yes | 0.08 | 0.00, 44.4 | 0.4 |
| Multi-level Injury | |||
| No | — | — | |
| Yes | 1,791 | 0.12, 175,968,272 | 0.14 |
| Side of Injury | |||
| Bilateral | — | — | |
| Left | 1.24 | 0.26, 6.77 | 0.8 |
| Right | 0.43 | 0.08, 2.51 | 0.3 |
| Mechanism of Injury | |||
| Sexual intercourse trauma | — | — | |
| Other | 0.65 | 0.00, 10.6 | 0.8 |
| Non-sexual blunt force trauma | 3.73 | 0.71, 18.2 | 0.12 |
| Masturbation injury | 5.43 | 0.63, 70.1 | 0.12 |
| Abbreviation: CI = Confidence Interval | |||
| Number of patients in analysis: 118 | |||
# 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_th5.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 < casesauc(pain_lr_roc)## Area under the curve: 0.8238ci.auc(pain_lr_roc)## 95% CI: 0.7176-0.9299 (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.104 |
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.346 |
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,
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 | |||
| Yes | — | — | |
| No | 2.65 | 1.02, 7.22 | 0.045 |
| Imaging | |||
| No | — | — | |
| Yes | 0.79 | 0.24, 2.58 | 0.7 |
| Centre Volume | |||
| High | — | — | |
| Medium | 1.65 | 0.37, 7.58 | 0.5 |
| Low | 1.51 | 0.42, 5.95 | 0.5 |
| Age at Time of Injury | 1.01 | 0.97, 1.06 | 0.6 |
| Urethral Injury | |||
| No | — | — | |
| Yes | 1.80 | 0.54, 5.98 | 0.3 |
| Transfer for Surgical Repair | |||
| No | — | — | |
| Yes | 1.52 | 0.46, 5.08 | 0.5 |
| Distal Corpus Location | |||
| No | — | — | |
| Yes | 31.2 | 0.06, 80,547 | 0.3 |
| Proximal Corpus Location | |||
| No | — | — | |
| Yes | 23.8 | 0.06, 53,304 | 0.3 |
| Mid Corpus Location | |||
| No | — | — | |
| Yes | 11.0 | 0.03, 24,877 | 0.5 |
| Crural location | |||
| No | — | — | |
| Yes | 6.29 | 0.03, 6,308 | 0.5 |
| Multi-level Injury | |||
| No | — | — | |
| Yes | 0.11 | 0.00, 192 | 0.6 |
| Side of Injury | |||
| Bilateral | — | — | |
| Left | 4.22 | 0.75, 46.5 | 0.11 |
| Right | 3.94 | 0.71, 43.2 | 0.12 |
| Mechanism of Injury | |||
| Sexual intercourse trauma | — | — | |
| Other | 0.67 | 0.00, 13.3 | 0.8 |
| Non-sexual blunt force trauma | 0.34 | 0.03, 1.74 | 0.2 |
| Masturbation injury | 0.22 | 0.00, 2.00 | 0.2 |
| Abbreviation: CI = Confidence Interval | |||
| Number of patients in analysis: 114 | |||
# 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 < casesauc(ed_lr_roc)## Area under the curve: 0.7761ci.auc(ed_lr_roc)## 95% CI: 0.6837-0.8684 (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.227 |
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,
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 +
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 = "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 |
|---|---|---|---|
| Theatre within 24 hours of Presentation | |||
| Yes | — | — | |
| No | 2.12 | 0.58, 8.39 | 0.3 |
| Imaging | |||
| No | — | — | |
| Yes | 0.33 | 0.06, 1.40 | 0.14 |
| Centre Volume | |||
| High | — | — | |
| Medium | 1.55 | 0.27, 10.3 | 0.6 |
| Low | 0.61 | 0.09, 4.30 | 0.6 |
| Age at Time of Injury | 1.00 | 0.94, 1.06 | >0.9 |
| Urethral Injury | |||
| No | — | — | |
| Yes | 1.56 | 0.30, 6.67 | 0.6 |
| Transfer for Surgical Repair | |||
| No | — | — | |
| Yes | 1.10 | 0.21, 5.13 | >0.9 |
| Site of Injury | |||
| Crura | — | — | |
| Distal | 0.70 | 0.01, 168 | 0.9 |
| Mid | 1.26 | 0.03, 281 | >0.9 |
| Multi-level | 1.07 | 0.00, 370 | >0.9 |
| Proximal | 1.33 | 0.03, 289 | 0.9 |
| Side of Injury | |||
| Bilateral | — | — | |
| Left | 0.25 | 0.03, 1.63 | 0.15 |
| Right | 0.76 | 0.14, 4.39 | 0.7 |
| Mechanism of Injury | |||
| Sexual intercourse trauma | — | — | |
| Other | 0.95 | 0.01, 23.0 | >0.9 |
| Non-sexual blunt force trauma | 0.20 | 0.00, 2.24 | 0.2 |
| Masturbation injury | 0.23 | 0.00, 3.35 | 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 < casesauc(waisting_shortening_lr_roc)## Area under the curve: 0.8286ci.auc(waisting_shortening_lr_roc)## 95% CI: 0.7409-0.9163 (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.102 |
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.209 |
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,
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 +
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 = "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 |
|---|---|---|---|
| Theatre within 24 hours of Presentation | |||
| Yes | — | — | |
| No | 1.71 | 0.57, 5.50 | 0.3 |
| Imaging | |||
| No | — | — | |
| Yes | 0.34 | 0.09, 1.17 | 0.087 |
| Centre Volume | |||
| High | — | — | |
| Medium | 0.57 | 0.11, 2.87 | 0.5 |
| Low | 0.41 | 0.09, 1.79 | 0.2 |
| Age at Time of Injury | 0.99 | 0.94, 1.04 | 0.7 |
| Urethral Injury | |||
| No | — | — | |
| Yes | 1.53 | 0.45, 5.08 | 0.5 |
| Transfer for Surgical Repair | |||
| No | — | — | |
| Yes | 3.14 | 0.94, 11.5 | 0.063 |
| Distal Corpus Location | |||
| No | — | — | |
| Yes | 0.34 | 0.00, 5.28 | 0.5 |
| Proximal Corpus Location | |||
| No | — | — | |
| Yes | 0.44 | 0.00, 9.51 | 0.6 |
| Mid Corpus Location | |||
| No | — | — | |
| Yes | 0.24 | 0.00, 4.23 | 0.4 |
| Crural location | |||
| No | — | — | |
| Yes | 0.29 | 0.00, 38.3 | 0.6 |
| laterality | |||
| Bilateral | — | — | |
| Left | 1.29 | 0.30, 5.67 | 0.7 |
| Right | 0.29 | 0.06, 1.22 | 0.090 |
| Mechanism of Injury | |||
| Sexual intercourse trauma | — | — | |
| Other | 0.38 | 0.00, 6.15 | 0.5 |
| Non-sexual blunt force trauma | 0.35 | 0.06, 1.72 | 0.2 |
| Masturbation injury | 0.15 | 0.00, 1.83 | 0.2 |
| 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 < casesauc(abnormal_curvature_lr_roc)## Area under the curve: 0.8449ci.auc(abnormal_curvature_lr_roc)## 95% CI: 0.7588-0.931 (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.142 |
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.38 |
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,
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 +
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 = "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.86 | 0.67, 1.04 | 0.13 |
| Imaging | |||
| No | — | — | |
| Yes | 1.28 | 0.39, 4.45 | 0.7 |
| Age at Time of Injury | 1.04 | 0.98, 1.10 | 0.2 |
| Urethral Injury | |||
| No | — | — | |
| Yes | 0.94 | 0.20, 3.58 | >0.9 |
| Centre Volume | |||
| High | — | — | |
| Medium | 8.46 | 1.20, 98.8 | 0.032 |
| Low | 11.7 | 2.09, 126 | 0.004 |
| Transfer for Surgical Repair | |||
| No | — | — | |
| Yes | 5.32 | 1.53, 25.6 | 0.008 |
| Distal Corpus Location | |||
| No | — | — | |
| Yes | 0.01 | 0.00, 65.5 | 0.4 |
| Proximal Corpus Location | |||
| No | — | — | |
| Yes | 0.01 | 0.00, 26.7 | 0.3 |
| Mid Corpus Location | |||
| No | — | — | |
| Yes | 0.01 | 0.00, 17.3 | 0.2 |
| Crural location | |||
| No | — | — | |
| Yes | 0.21 | 0.00, 121 | 0.6 |
| Multi-level Injury | |||
| No | — | — | |
| Yes | 505 | 0.04, 24,621,818 | 0.2 |
| Side of Injury | |||
| Bilateral | — | — | |
| Left | 0.55 | 0.11, 2.81 | 0.5 |
| Right | 0.25 | 0.05, 1.35 | 0.10 |
| Mechanism of Injury | |||
| Sexual intercourse trauma | — | — | |
| Other | 0.50 | 0.00, 8.25 | 0.7 |
| Non-sexual blunt force trauma | 3.84 | 0.73, 19.1 | 0.11 |
| Masturbation injury | 8.59 | 0.88, 119 | 0.065 |
| Abbreviation: CI = Confidence Interval | |||
| Number of patients in analysis: 119 | |||
# 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_th6.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 < casesauc(pain_lr_roc)## Area under the curve: 0.5504ci.auc(pain_lr_roc)## 95% CI: 0.4312-0.6696 (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.104 |
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.352 |
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 +
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 = "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.39 | 0.11 |
| Imaging | |||
| No | — | — | |
| Yes | 0.67 | 0.19, 2.34 | 0.5 |
| Centre Volume | |||
| High | — | — | |
| Medium | 1.61 | 0.34, 7.79 | 0.5 |
| Low | 1.36 | 0.36, 5.67 | 0.7 |
| Age at Time of Injury | 1.02 | 0.97, 1.07 | 0.5 |
| Urethral Injury | |||
| No | — | — | |
| Yes | 2.55 | 0.73, 9.24 | 0.14 |
| Transfer for Surgical Repair | |||
| No | — | — | |
| Yes | 1.42 | 0.43, 4.77 | 0.6 |
| Distal Corpus Location | |||
| No | — | — | |
| Yes | 9.20 | 0.01, 30,261 | 0.5 |
| Proximal Corpus Location | |||
| No | — | — | |
| Yes | 13.1 | 0.02, 33,764 | 0.4 |
| Mid Corpus Location | |||
| No | — | — | |
| Yes | 4.32 | 0.01, 11,483 | 0.7 |
| Crural location | |||
| No | — | — | |
| Yes | 6.33 | 0.03, 6,959 | 0.5 |
| Multi-level Injury | |||
| No | — | — | |
| Yes | 0.27 | 0.00, 685 | 0.8 |
| Side of Injury | |||
| Bilateral | — | — | |
| Left | 6.78 | 1.16, 82.6 | 0.032 |
| Right | 4.43 | 0.73, 55.2 | 0.11 |
| Mechanism of Injury | |||
| Sexual intercourse trauma | — | — | |
| Other | 0.42 | 0.00, 8.92 | 0.6 |
| Non-sexual blunt force trauma | 0.35 | 0.04, 1.78 | 0.2 |
| Masturbation injury | 0.38 | 0.00, 4.28 | 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 < casesauc(ed_lr_roc)## Area under the curve: 0.5014ci.auc(ed_lr_roc)## 95% CI: 0.4216-0.5812 (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 lengthcbind(
"Score" = "Brier Score",
"Brier Score" = round(brier_score, digits = 3)
) %>% as_tibble() %>% gt() %>% gt_theme_espn()| Score | Brier Score |
|---|---|
| Brier Score | 0.213 |
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.225 |
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,
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 +
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 = "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.94 | 0.65, 1.28 | 0.7 |
| Imaging | |||
| No | — | — | |
| Yes | 0.31 | 0.05, 1.59 | 0.2 |
| Centre Volume | |||
| High | — | — | |
| Medium | 0.92 | 0.12, 7.25 | >0.9 |
| Low | 0.36 | 0.03, 3.18 | 0.4 |
| Age at Time of Injury | 1.03 | 0.95, 1.13 | 0.4 |
| Urethral Injury | |||
| No | — | — | |
| Yes | 2.88 | 0.53, 14.3 | 0.2 |
| Transfer for Surgical Repair | |||
| No | — | — | |
| Yes | 0.95 | 0.17, 4.55 | >0.9 |
| Site of Injury | |||
| Crura | — | — | |
| Distal | 0.04 | 0.00, 17.9 | 0.3 |
| Mid | 0.19 | 0.00, 45.6 | 0.5 |
| Multi-level | 0.32 | 0.00, 102 | 0.7 |
| Proximal | 0.32 | 0.00, 75.9 | 0.6 |
| Side of Injury | |||
| Bilateral | — | — | |
| Left | 0.25 | 0.03, 1.62 | 0.14 |
| Right | 0.44 | 0.07, 2.75 | 0.4 |
| Mechanism of Injury | |||
| Sexual intercourse trauma | — | — | |
| Other | 0.48 | 0.00, 12.0 | 0.7 |
| Non-sexual blunt force trauma | 0.42 | 0.00, 5.04 | 0.5 |
| Masturbation injury | 3.38 | 0.01, 144 | 0.6 |
| Abbreviation: CI = Confidence Interval | |||
| Number of patients in analysis: 98 | |||
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 > casesauc(waisting_shortening_lr_roc)## Area under the curve: 0.5689ci.auc(waisting_shortening_lr_roc)## 95% CI: 0.4015-0.7364 (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.084 |
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.254 |
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,
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 +
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 = "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.02 | 0.84, 1.24 | 0.8 |
| Imaging | |||
| No | — | — | |
| Yes | 0.33 | 0.08, 1.21 | 0.10 |
| Centre Volume | |||
| High | — | — | |
| Medium | 0.71 | 0.13, 3.87 | 0.7 |
| Low | 0.38 | 0.08, 1.81 | 0.2 |
| Age at Time of Injury | 0.99 | 0.94, 1.05 | 0.8 |
| Urethral Injury | |||
| No | — | — | |
| Yes | 1.48 | 0.39, 5.31 | 0.5 |
| Transfer for Surgical Repair | |||
| No | — | — | |
| Yes | 3.08 | 0.91, 11.9 | 0.072 |
| Distal Corpus Location | |||
| No | — | — | |
| Yes | 0.40 | 0.00, 6.62 | 0.6 |
| Proximal Corpus Location | |||
| No | — | — | |
| Yes | 0.42 | 0.00, 9.90 | 0.6 |
| Mid Corpus Location | |||
| No | — | — | |
| Yes | 0.18 | 0.00, 3.32 | 0.3 |
| Crural location | |||
| No | — | — | |
| Yes | 0.43 | 0.00, 68.7 | 0.8 |
| Side of Injury | |||
| Bilateral | — | — | |
| Left | 1.28 | 0.29, 5.86 | 0.7 |
| Right | 0.26 | 0.05, 1.16 | 0.078 |
| Mechanism of Injury | |||
| Sexual intercourse trauma | — | — | |
| Other | 0.31 | 0.00, 4.61 | 0.4 |
| Non-sexual blunt force trauma | 0.32 | 0.05, 1.69 | 0.2 |
| Masturbation injury | 0.37 | 0.00, 4.43 | 0.5 |
| Abbreviation: CI = Confidence Interval | |||
| Number of patients in analysis: 102 | |||
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 < casesauc(abnormal_curvature_lr_roc)## Area under the curve: 0.5753ci.auc(abnormal_curvature_lr_roc)## 95% CI: 0.4697-0.6809 (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.143 |
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.407 |