library(tidyverse)
library(tidymodels)
library(magrittr)
library(here) # for path management
library(readr)
library(openxlsx)
library(readr) # for importing data
library(broom.mixed) # for converting Bayesian models to tidy tibbles
library(dotwhisker) # for visualizing regression results
library(skimr) # for variable summaries
library(tableone) # for tables
library(gtsummary) # for prettier tables
library(patchwork) # for multi-planel figures
library(RColorBrewer) # for color palettesJones Fractures
data_path <- here('5mt_data_20230507.xlsx')
df.raw <- read.xlsx(data_path,
detectDates = TRUE)
df.raw %>% colnames() [1] "physician.1" "physician.2" "specialty"
[4] "name" "notes" "side"
[7] "mrn" "zone" "jones"
[10] "base" "include" "sex"
[13] "race" "smoker" "bmi"
[16] "age" "comorb" "cci"
[19] "diabetes" "ambulation_baseline" "device_baseline"
[22] "date_xr_initial" "date_presentation" "date_treating_office"
[25] "date_final" "no_postacute_fu" "mechanism"
[28] "injuries" "initial_pres_uc" "uc_wb"
[31] "uc_immob" "initial_pres_ed" "ed_wb"
[34] "ed_immob" "initial_pres_office" "office_wb"
[37] "office_immob" "operative" "fixation"
[40] "implants" "failed_nonop" "date_failed_nonop"
[43] "splint" "date_splint_on" "date_splint_off"
[46] "cam" "date_cam_on" "date_cam_off"
[49] "hardsole" "date_hardsole_on" "date_hardsole_off"
[52] "date_rts" "nonunion" "refracture"
[55] "infection" "nerve_inj" "device_at_final_fu"
[58] "comments" "date_last_nonhealed_xr" "date_injury"
[61] "date_or" "date_clinical_healed" "date_rad_healed" # Encode qualitative columns as factors
df.raw <- df.raw %>% mutate_if(is.character, as.factor)
# Fix variable types
df.raw <- df.raw %>%
mutate(
across(c(mrn,
zone,
cci),
factor),
across(starts_with('date_'),
ymd)
)
# New vars
df.raw <- df.raw %>%
mutate(duration_rad_healed = case_when(
failed_nonop == 'No' ~ difftime(date_rad_healed, date_injury, units = 'weeks') %>% as.numeric(),
is.na(failed_nonop) ~ difftime(date_rad_healed, date_injury, units = 'weeks') %>% as.numeric(),
failed_nonop == 'Yes' ~ difftime(date_rad_healed, date_or, units = 'weeks') %>% as.numeric()
)) %>%
mutate(duration_clinical_healed = case_when(
failed_nonop == 'No' ~ difftime(date_clinical_healed, date_injury, units = 'weeks') %>% as.numeric(),
is.na(failed_nonop) ~ difftime(date_clinical_healed, date_injury, units = 'weeks') %>% as.numeric(),
failed_nonop == 'Yes' ~ difftime(date_clinical_healed, date_or, units = 'weeks') %>% as.numeric()
)) %>%
mutate(duration_injury_to_office = difftime(date_treating_office, date_injury, units = 'weeks') %>% as.numeric()) %>%
mutate(duration_trial_of_nonop = difftime(date_or, date_injury, units = 'weeks') %>% as.numeric()) %>%
mutate(duration_rts = difftime(date_rts, date_injury, units = 'weeks') %>% as.numeric()) %>%
mutate(followup = difftime(date_final, date_presentation, units = 'weeks') %>% as.numeric()) %>%
mutate(duration_immob_splint = difftime(date_splint_off, date_splint_on, units = 'weeks') %>% as.numeric()) %>%
mutate(duration_immob_cam = difftime(date_cam_off, date_cam_on, units = 'weeks') %>% as.numeric()) %>%
mutate(duration_immob_hardsole = difftime(date_hardsole_off, date_hardsole_on, units = 'weeks') %>% as.numeric()) %>%
mutate(duration_immobilization = select(., all_of(starts_with('duration_immob_'))) %>% rowSums(., na.rm = TRUE)) %>%
mutate(duration_immobilization = case_when(
duration_immobilization > 0 ~ duration_immobilization,
duration_immobilization == 0 ~ NA
)) %>%
mutate(initial_presentation_setting = case_when(
initial_pres_ed == 'Yes' ~ 'ed',
initial_pres_uc == 'Yes' ~ 'uc',
initial_pres_office == 'Yes' ~ 'office',
TRUE ~ NA_character_
)) %>%
mutate(initial_wb = case_when(
initial_presentation_setting == 'ed' ~ ed_wb,
initial_presentation_setting == 'uc' ~ uc_wb,
initial_presentation_setting == 'office' ~ office_wb,
TRUE ~ NA_character_
)) %>%
mutate(splint = case_when(
splint == 'Yes' ~ 'Yes',
is.na(splint) ~ 'No'
)) %>%
mutate(cam = case_when(
cam == 'Yes' ~ 'Yes',
is.na(cam) ~ 'No'
)) %>%
mutate(hardsole = case_when(
hardsole == 'Yes' ~ 'Yes',
is.na(hardsole) ~ 'No'
)) %>%
mutate(failed_nonop = case_when(
failed_nonop == 'Yes' ~ 'Yes',
is.na(failed_nonop) ~ 'No'
)) %>%
mutate(duration_without_healing = case_when(
!is.na(date_rad_healed) ~ NA,
is.na(date_rad_healed) ~ difftime(date_last_nonhealed_xr, date_injury, units = 'weeks') %>% as.numeric()
)) %>%
mutate(delayed_union = case_when(
duration_rad_healed*7 > 180 ~ 'Yes',
duration_rad_healed*7 <= 180 ~ 'No',
is.na(duration_rad_healed) & duration_without_healing*7 <= 180 ~ 'No',
is.na(duration_rad_healed) & duration_without_healing*7 > 180 ~ 'Yes'
)) %>%
mutate(failed_nonop_plus_delayed_union = case_when(
operative == 'Nonoperative' & delayed_union == 'Yes' ~ 'Yes',
failed_nonop == 'Yes' ~ 'Yes',
TRUE ~ 'No'
))
# Create initial treatment var
df.raw <- df.raw %>%
mutate(initial_treatment = case_when(
operative == 'Nonoperative' ~ 'nonop',
operative == 'Operative' & failed_nonop == 'Yes' ~ 'nonop',
operative == 'Operative' & failed_nonop == 'No' ~ 'operative',
TRUE ~ NA_character_
)) %>%
mutate(initial_treatment = initial_treatment %>% as.factor())
# Exclude patients
df.raw <- df.raw %>% filter(include != '0')
# # Reorder and rename factors/levels
# df.raw <- df.raw %>%
# mutate(var1 =
# factor(var1,
# levels = c('level1',
# 'level2',
# 'level3')))
df <- df.raw %>%
filter(jones == 'Y')
# df of jones fractures
df <- df %>%
filter(jones == 'Y',
!is.na(operative))n_patients_all <- df.raw %>% distinct(mrn) %>% nrow() %>% as.numeric()
n_jones_all <- df.raw %>% filter(jones == 'Y') %>% nrow() %>% as.numeric()
n_jones <- df %>% filter(!is.na(operative)) %>% nrow() %>% as.numeric()
n_initial_nonoperative <- df %>% filter(initial_treatment == 'nonop') %>% nrow() %>% as.numeric()
n_initial_operative <- df %>% filter(initial_treatment == 'operative') %>% nrow() %>% as.numeric()
# Demographics
demos <- df %>%
filter(!is.na(operative)) %>%
summarise(
across(
c(age,
bmi),
list(mean = \(x) mean(x, na.rm = TRUE),
sd = \(x) sd(x, na.rm = TRUE),
min = \(x) min(x, na.rm = TRUE),
max = \(x) max(x, na.rm = TRUE)),
.names = '{.col}.{.fn}'
)
) %>%
pivot_longer(everything()) %>%
mutate(
value = round(value, digits = 1)
) %>%
separate_wider_delim(
cols = name,
delim = '.',
names = c('var', 'stat')
) %>%
pivot_wider(
names_from = stat,
values_from = value
) %>%
as.data.frame() %>%
column_to_rownames(var = 'var')
# # Followup
# demos <- df %>%
# mutate(fu_fix = case_when(
# followup < 1.7 & !is.na(MHHS_2yr) ~ 2,
# TRUE ~ followup
# )) %>%
# filter(case_id == '1') %>%
# filter(followup >= 1.7) %>%
# summarise(
# across(c(fu_fix),
# list(mean = mean,
# sd = sd),
# .names = '{.col}.{.fn}')
# ) %>%
# pivot_longer(everything()) %>%
# mutate(
# value = round(value, digits = 1)
# ) %>%
# separate_wider_delim(
# cols = name,
# delim = '.',
# names = c('var', 'stat')
# ) %>%
# pivot_wider(
# names_from = stat,
# values_from = value
# ) %>%
# as.data.frame() %>%
# column_to_rownames(var = 'var') %>%
# bind_rows(demos)
# Nonop failed
n_failed_nonop <- df %>% filter(failed_nonop == 'Yes') %>% nrow() %>% as.numeric()
demos_failed_nonop <- df %>%
filter(failed_nonop == 'Yes') %>%
summarise(
across(
c(duration_trial_of_nonop),
list(mean = mean,
sd = sd,
min = min,
max = max),
.names = '{.col}.{.fn}'
)
) %>%
pivot_longer(everything()) %>%
mutate(
value = round(value, digits = 1)
) %>%
separate_wider_delim(
cols = name,
delim = '.',
names = c('var', 'stat')
) %>%
pivot_wider(
names_from = stat,
values_from = value
) %>%
as.data.frame() %>%
column_to_rownames(var = 'var')Results
Demographics
A total of 2450 patients presented with a 5th metatarsal fracture during the study period. Among all 5th metatarsal fractures, a total of 165 fractures (6.7%) were diagnosed as a true Jones fracture based on radiographs. 45 patients did not follow up within our health system after the initial visit and were excluded from this analysis. Therefore, a total of 120 patients with Jones fractures were included in this study. The mean age of the overall cohort was 46.7 +/- 18.5 years and the mean BMI was 28.2 +/- 6.2.
Among the included patients, 99 patients (82.5%) were initially treated nonoperatively and 21 patients (17.5%) were treated operatively. While there were differences in patient race between groups, there were otherwise no differences in baseline demographic factors between patients treated operatively and those treated nonoperatively (Table 1).
t1 <- df %>%
filter(!is.na(operative)) %>%
select(
age,
sex,
bmi,
race,
smoker,
cci,
followup,
initial_treatment
) %>%
mutate(initial_treatment = recode_factor(initial_treatment,
'nonop' = 'Nonoperative',
'operative' = 'Operative'
)) %>%
droplevels %>%
tbl_summary(
by = initial_treatment,
missing = 'no',
statistic = list(
all_continuous() ~ '{mean} +/- {sd}'
),
digits = list(
all_categorical() ~ c(0,1),
all_continuous() ~ c(1,1)),
label = list(
age ~ 'Age (years)',
sex ~ 'Sex',
bmi ~ 'BMI',
race ~ 'Race',
smoker ~ 'Smoking status',
cci ~ 'Charlson Comorbidity Index',
followup ~ 'Follow-up duration (weeks)'
)
) %>%
add_p(test = list(
all_continuous() ~ 't.test'
),
pvalue_fun = function(x) style_pvalue(x, digits = 3),
test.args = all_tests('fisher.test') ~ list(simulate.p.value = TRUE)) %>%
bold_p(t = 0.05) %>%
add_overall() %>%
modify_spanning_header(all_stat_cols(stat_0 = F) ~ "**Initial Treatment**") %>%
modify_caption("<div style='text-align: left; font-weight: bold'>
Overall Patient Demographics") %>%
bold_labels() %>%
as_gt() %>%
gt::tab_options(
table.font.size = 'small',
data_row.padding = gt::px(3)
)
t1| Characteristic | Overall, N = 1201 | Initial Treatment | p-value2 | |
|---|---|---|---|---|
| Nonoperative, N = 991 | Operative, N = 211 | |||
| Age (years) | 46.7 +/- 18.5 | 47.4 +/- 18.7 | 43.3 +/- 17.5 | 0.345 |
| Sex | 0.225 | |||
| Female | 82 (68.3%) | 70 (70.7%) | 12 (57.1%) | |
| Male | 38 (31.7%) | 29 (29.3%) | 9 (42.9%) | |
| BMI | 28.2 +/- 6.2 | 28.2 +/- 6.1 | 28.4 +/- 7.0 | 0.891 |
| Race | 0.006 | |||
| Asian | 11 (9.2%) | 11 (11.1%) | 0 (0.0%) | |
| Black | 13 (10.8%) | 9 (9.1%) | 4 (19.0%) | |
| Other | 17 (14.2%) | 9 (9.1%) | 8 (38.1%) | |
| Pacific Islander | 1 (0.8%) | 1 (1.0%) | 0 (0.0%) | |
| Unknown | 14 (11.7%) | 14 (14.1%) | 0 (0.0%) | |
| White | 64 (53.3%) | 55 (55.6%) | 9 (42.9%) | |
| Smoking status | 0.736 | |||
| Daily | 13 (10.9%) | 12 (12.2%) | 1 (4.8%) | |
| Former | 19 (16.0%) | 16 (16.3%) | 3 (14.3%) | |
| Never | 87 (73.1%) | 70 (71.4%) | 17 (81.0%) | |
| Charlson Comorbidity Index | 0.556 | |||
| 0 | 91 (76.5%) | 76 (77.6%) | 15 (71.4%) | |
| 1 | 17 (14.3%) | 14 (14.3%) | 3 (14.3%) | |
| 2 | 8 (6.7%) | 5 (5.1%) | 3 (14.3%) | |
| 3 | 2 (1.7%) | 2 (2.0%) | 0 (0.0%) | |
| 6 | 1 (0.8%) | 1 (1.0%) | 0 (0.0%) | |
| Follow-up duration (weeks) | 18.3 +/- 17.9 | 17.0 +/- 16.8 | 24.4 +/- 21.8 | 0.155 |
| 1 Mean +/- SD; n (%) | ||||
| 2 Welch Two Sample t-test; Pearson’s Chi-squared test; Fisher’s Exact Test for Count Data with simulated p-value (based on 2000 replicates) | ||||
Patients treated nonoperatively
df_nonop <- df %>%
filter(initial_treatment == 'nonop')
# perc_nonop_ed <- df_nonop %>%
# filter(initial_presentation_setting == 'ed') %>%
# nrow() %>%
# as.numeric() %>%
# divide_by(n_nonoperative) %>%
# multiply_by(100) %>%
# round(digits = 1)
# perc_nonop_uc <- df_nonop %>%
# filter(initial_presentation_setting == 'uc') %>%
# nrow() %>%
# as.numeric() %>%
# divide_by(n_nonoperative) %>%
# multiply_by(100) %>%
# round(digits = 1)
# perc_nonop_office <- df_nonop %>%
# filter(initial_presentation_setting == 'office') %>%
# nrow() %>%
# as.numeric() %>%
# divide_by(n_nonoperative) %>%
# multiply_by(100) %>%
# round(digits = 1)
nonop_summary <- df_nonop %>%
summarise(
across(
c(duration_immob_splint,
duration_immob_cam,
duration_immob_hardsole,
duration_immobilization),
list(mean = \(x) mean(x, na.rm = TRUE),
sd = \(x) sd(x, na.rm = TRUE),
min = \(x) min(x, na.rm = TRUE),
max = \(x) max(x, na.rm = TRUE)),
.names = '{.col}.{.fn}'
)
) %>%
pivot_longer(everything()) %>%
mutate(
value = round(value, digits = 1)
) %>%
separate_wider_delim(
cols = name,
delim = '.',
names = c('var', 'stat')
) %>%
pivot_wider(
names_from = stat,
values_from = value
) %>%
as.data.frame() %>%
column_to_rownames(var = 'var')There was considerable variability in the treatment courses among patients treated without surgery (Table 2). 35.7% of patients presented initially to an emergency department, 36.7% presented to an urgent care center, and 27.6% presented directly to a treating physician’s office. Upon initial presentation, 74.7% of patients were made non-weightbearing, while 5.1% were made partial weightbearing and 20.3% were allowed to bear weight as tolerated. After being evaluated in the office, 61.4% of patients were made non-weightbearing, 9.1% partial weightbearing, and 29.5% weightbearing as tolerated. Over the duration of the treatment course, 36.4% of patients were placed in a splint or cast for a mean duration of 1.9 +/- 2 weeks. 81.8% of patients were treated with a controlled ankle motion (CAM) boot for a mean duration of 8.2 +/- 4.5 weeks, either acutely upon initial presentation or following splint/cast immobilization. Likewise, 33.3% of patients were further managed with a hard sole shoe during their treatment course for a mean duration of 4.2 +/- 4 weeks before transitioning to normal footwear. Overall, patients treated nonoperatively were managed with some form of foot immobilization or support (cast, splint, CAM boot, or hard sole shoe) for an average of 9 +/- 4.8 weeks.
4 patients were initially treated nonoperatively but later underwent operative management for delayed healing (4% of patients treated nonoperatively). The mean duration of failed nonoperative treatment before undergoing surgery was 25 +/- 21 weeks (range 8 to 53.9 weeks).
t2 <- df %>%
filter(initial_treatment == 'nonop') %>%
select(
initial_presentation_setting,
initial_wb,
office_wb,
splint,
cam,
hardsole,
duration_immobilization
) %>%
mutate(initial_presentation_setting = recode_factor(initial_presentation_setting,
'ed' = 'Emergency department',
'uc' = 'Urgent care',
'office' = 'Office'
)) %>%
mutate(initial_wb = recode_factor(initial_wb,
'NWB' = 'Non-weightbearing',
'PWB' = 'Partial weightbearing',
'WBAT' = 'Weightbearing as tolerated'
)) %>%
mutate(office_wb = recode_factor(office_wb,
'NWB' = 'Non-weightbearing',
'PWB' = 'Partial weightbearing',
'WBAT' = 'Weightbearing as tolerated'
)) %>%
droplevels %>%
tbl_summary(
missing = 'no',
statistic = list(
all_continuous() ~ '{mean} +/- {sd}',
all_categorical() ~ '{p}%'
),
digits = list(
all_categorical() ~ c(1),
all_continuous() ~ c(1,1)),
label = list(
initial_presentation_setting ~ 'Initial presentation setting',
initial_wb ~ 'Initial weightbearing restriction',
office_wb ~ 'Treating office weightbearing restriction',
splint ~ 'Splint or cast',
cam ~ 'CAM or walking boot',
hardsole ~ 'Hard sole shoe',
duration_immobilization ~ 'Duration of foot immobilization (weeks)'
)
) %>%
modify_caption("<div style='text-align: left; font-weight: bold'>
Treatment course for patients initially treated nonoperatively") %>%
bold_labels() %>%
as_gt() %>%
gt::tab_options(
table.font.size = 'small',
data_row.padding = gt::px(3)
)
t2| Characteristic | N = 991 |
|---|---|
| Initial presentation setting | |
| Emergency department | 35.7% |
| Urgent care | 36.7% |
| Office | 27.6% |
| Initial weightbearing restriction | |
| Non-weightbearing | 74.7% |
| Partial weightbearing | 5.1% |
| Weightbearing as tolerated | 20.3% |
| Treating office weightbearing restriction | |
| Non-weightbearing | 61.4% |
| Partial weightbearing | 9.1% |
| Weightbearing as tolerated | 29.5% |
| Splint or cast | 36.4% |
| CAM or walking boot | 81.8% |
| Hard sole shoe | 33.3% |
| Duration of foot immobilization (weeks) | 9.0 +/- 4.8 |
| 1 %; Mean +/- SD | |
Outcomes
outcomes_summary <- df %>%
group_by(operative) %>%
summarise(
across(
c(duration_immobilization,
duration_rad_healed,
duration_clinical_healed),
list(mean = \(x) mean(x, na.rm = TRUE),
sd = \(x) sd(x, na.rm = TRUE),
min = \(x) min(x, na.rm = TRUE),
max = \(x) max(x, na.rm = TRUE)),
.names = '{.col}.{.fn}'
)
) %>%
pivot_longer(-operative) %>%
mutate(
value = round(value, digits = 1)
) %>%
separate_wider_delim(
cols = name,
delim = '.',
names = c('var', 'stat')
) %>%
pivot_wider(
id_cols = c(operative, var),
names_from = stat,
values_from = value,
) %>%
unite(
'var',
operative:var,
sep = '_'
) %>%
as.data.frame() %>%
column_to_rownames(var = 'var')
t.test(duration_immobilization ~ operative, data = df)
Welch Two Sample t-test
data: duration_immobilization by operative
t = 0.6531, df = 29.138, p-value = 0.5188
alternative hypothesis: true difference in means between group Nonoperative and group Operative is not equal to 0
95 percent confidence interval:
-1.827899 3.543500
sample estimates:
mean in group Nonoperative mean in group Operative
8.702521 7.844720 pval1 <- 0.519
t.test(duration_clinical_healed ~ operative, data = df)
Welch Two Sample t-test
data: duration_clinical_healed by operative
t = 0.25396, df = 47.312, p-value = 0.8006
alternative hypothesis: true difference in means between group Nonoperative and group Operative is not equal to 0
95 percent confidence interval:
-2.281556 2.940954
sample estimates:
mean in group Nonoperative mean in group Operative
12.66541 12.33571 pval2 <- 0.801
t.test(duration_rad_healed ~ operative, data = df)
Welch Two Sample t-test
data: duration_rad_healed by operative
t = 1.0385, df = 46.879, p-value = 0.3044
alternative hypothesis: true difference in means between group Nonoperative and group Operative is not equal to 0
95 percent confidence interval:
-1.563636 4.900201
sample estimates:
mean in group Nonoperative mean in group Operative
13.20724 11.53896 pval3 <- 0.304There was no evidence for any difference between groups with respect to the time until clinical healing (12.7 +/- 7.1 weeks for the nonoperative group versus 12.3 +/- 4.5 weeks for the operative group, p = 0.801; Table 3 and Figure 1). Likewise, there was no difference between groups with respect to the time until radiographic healing (13.2 +/- 8.1 weeks for the nonoperative group versus 11.5 +/- 6 weeks for the operative group, p = 0.304). In the nonoperative group, 4 patients (4.2%) experienced delayed radiographic union with time to healing greater than 180 days. In the operative group, 1 patient (4.0%) experienced delayed radiographic union.
There was no difference between groups with respect to the total amount of time spent with the foot immobilized or supported before transitioning to normal footwear (8.7 +/- 4.4 weeks for the nonoperative group versus 7.8 +/- 5.9 weeks for the operative group, p = 0.519).
t3 <- df %>%
mutate(delayed_union = case_when(
!is.na(delayed_union) ~ delayed_union,
is.na(delayed_union) ~ 'No'
)) %>%
select(
duration_clinical_healed,
duration_rad_healed,
delayed_union,
operative
) %>%
# mutate(initial_treatment = recode_factor(initial_treatment,
# 'nonop' = 'Nonoperative',
# 'operative' = 'Operative'
# )) %>%
droplevels %>%
tbl_summary(
by = operative,
missing = 'no',
statistic = list(
all_continuous() ~ '{mean} +/- {sd}'
),
digits = list(
all_categorical() ~ c(0,1),
all_continuous() ~ c(1,1)),
label = list(
duration_clinical_healed ~ 'Time to clinical healing (weeks)',
duration_rad_healed ~ 'Time to radiographic healing (weeks)',
delayed_union ~ 'Delayed union (> 180 days)'
)
) %>%
add_p(test = list(
all_continuous() ~ 't.test'
),
pvalue_fun = function(x) style_pvalue(x, digits = 3),
test.args = all_tests('fisher.test') ~ list(simulate.p.value = TRUE)) %>%
bold_p(t = 0.05) %>%
modify_spanning_header(all_stat_cols(stat_0 = F) ~ "**Treatment**") %>%
modify_caption("<div style='text-align: left; font-weight: bold'>
Outcomes") %>%
bold_labels() %>%
as_gt() %>%
gt::tab_options(
table.font.size = 'small',
data_row.padding = gt::px(3)
)
t3| Characteristic | Treatment | p-value2 | |
|---|---|---|---|
| Nonoperative, N = 951 | Operative, N = 251 | ||
| Time to clinical healing (weeks) | 12.7 +/- 7.1 | 12.3 +/- 4.5 | 0.801 |
| Time to radiographic healing (weeks) | 13.2 +/- 8.1 | 11.5 +/- 6.0 | 0.304 |
| Delayed union (> 180 days) | 4 (4.2%) | 1 (4.0%) | >0.999 |
| 1 Mean +/- SD; n (%) | |||
| 2 Welch Two Sample t-test; Fisher’s exact test | |||
library(ggridges)
# Radiographic healing
d_max <- density(df$duration_rad_healed, na.rm = TRUE)
d_min <- range(df$duration_rad_healed, na.rm = TRUE)
m <- df %>%
group_by(operative) %>%
summarise(median = median(duration_rad_healed, na.rm = TRUE))
p1 <- df %>%
ggplot(aes(x = duration_rad_healed,
fill = operative)) +
geom_density(alpha = 0.5,
size = 0.2) +
geom_vline(data = m,
aes(xintercept = median,
color = operative),
linetype = 'dashed',
size = 0.5) +
scale_x_continuous(breaks = seq(0,52,4),
limits = c(min(d_min), max(d_max$x))) +
scale_y_continuous(breaks = seq(0,0.12, 0.03),
limits = c(0,0.12)) +
labs(x = 'Weeks until radiographic healing',
y = 'Density') +
scale_fill_brewer(palette = 'Set1',
name = 'Treatment') +
scale_color_brewer(palette = 'Set1',
name = 'Treatment')
# Clinical healing
d_max2 <- density(df$duration_clinical_healed, na.rm = TRUE)
d_min2 <- range(df$duration_clinical_healed, na.rm = TRUE)
m2 <- df %>%
group_by(operative) %>%
summarise(median = median(duration_clinical_healed, na.rm = TRUE))
p2 <- df %>%
ggplot(aes(x = duration_clinical_healed,
fill = operative)) +
geom_density(alpha = 0.5,
size = 0.2) +
geom_vline(data = m2,
aes(xintercept = median,
color = operative),
linetype = 'dashed',
size = 0.5) +
scale_x_continuous(breaks = seq(0,52,4),
limits = c(min(d_min2), max(d_max$x))) +
scale_y_continuous(breaks = seq(0,0.12, 0.03),
limits = c(0,0.12)) +
labs(x = 'Weeks until clinical healing',
y = 'Density') +
scale_fill_brewer(palette = 'Set1',
name = 'Treatment') +
scale_color_brewer(palette = 'Set1',
name = 'Treatment')
fig1 <- (p1 / p2) + plot_annotation(tag_levels = 'A')
fig1
ggsave('5mt-jones-fractures_fig1.png',
path = 'C:/Users/kinge/Dropbox/Research/Projects/Jones fractures/5th MT Manuscripts/Jones op vs nonop',
plot = fig1,
dpi = 320,
width = 7.5,
height = 5.7,
units = 'in')