Effect of the use of 3D printed replicas on the surgical time of autotransplantation: statistical analysis

Published

November 9, 2022

Show the code

knitr::opts_chunk$set(warning = FALSE, message = FALSE)

If not installed, install pacman package

Show the code

pacman::p_load(tidyverse, 
               gtsummary, 
               janitor,
               sjPlot, 
               gt, 
               gtExtras, 
               tidymodels, 
               effects, 
               robustbase, 
               ggpval, # p values for figures 
               easystats)

Show the code

theme_set(theme_minimal())

Dataset

Show the code

df %>% 
  slice_sample(n = 6) %>% 
  knitr::kable()

patient_nr	gender	group	age_surgery	total_surgery_time_min	extra_alveolar_time_s	donor_fitting_times	stage_of_develop	donor_tooth_recipient_region	x10	date_of_birth	date_of_surgery	last_control	outcome	average_age	a
38	Female	1	14	35	30	1X	3	18-16	NA	1/24/2007	10/5/2021	9/22/2022	successful	NA	NA
44	Female	1	16	60	60	3x	4	28- 16	NA	5/3/2006	9/20/2022	9/22/2022	too early to say, but surgery was successful	NA	NA
16	Female	0	20	90	60	2X	4	38-46	NA	2/18/2000	7/28/2020	9/7/2022	successful	NA	NA
42	Male	1	17	40	60	2x	3	28-16	NA	9/28/2004	5/23/2022	6/16/2022	successful	NA	NA
1	Female	0	20	50	180	3x	3	28- 36	NA	11/26/1998	8/14/2019	8/14/2021	successful	NA	NA
13	Male	0	18	120	30	1X	3	48-37	NA	5/4/2002	6/29/2020	7/7/2021	extracted	PA infection	NA

Data cleaning

Show the code

df <- df %>% 
  select(-x10, 
         -a,
         patient_nr, 
         -average_age) %>% 
  mutate(donor_fitting_times = readr::parse_number(donor_fitting_times) , 
         date_of_birth = lubridate::mdy(date_of_birth), 
         date_of_surgery = lubridate::mdy(date_of_surgery), 
         last_control = lubridate::mdy(last_control), 
         donor_tooth_recipient_region = str_replace(donor_tooth_recipient_region, " ", ""))

Show the code

glimpse(df)

Rows: 46
Columns: 13
$ patient_nr                   <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13…
$ gender                       <chr> "Female", "Male", "Female", "Male", "Fema…
$ group                        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
$ age_surgery                  <dbl> 20, 19, 21, 17, 20, 19, 18, 16, 14, 17, 1…
$ total_surgery_time_min       <dbl> 50, 30, 60, 70, 110, 70, 50, 60, 50, 70, …
$ extra_alveolar_time_s        <dbl> 180, 30, 30, 90, 180, 30, 30, 120, 60, 30…
$ donor_fitting_times          <dbl> 3, 1, 2, 3, 5, 1, 3, 3, 2, 1, 2, 2, 1, 1,…
$ stage_of_develop             <dbl> 3, 5, 5, 4, 4, 5, 4, 4, 3, 3, 5, 3, 3, 4,…
$ donor_tooth_recipient_region <chr> "28-36", "38-37", "38-36", "18-36", "38-1…
$ date_of_birth                <date> 1998-11-26, 2000-08-01, 1998-06-17, 2001…
$ date_of_surgery              <date> 2019-08-14, 2019-08-27, 2019-12-09, 2019…
$ last_control                 <date> 2021-08-14, 2021-08-12, 2022-09-07, 2021…
$ outcome                      <chr> "successful", "successful", "successful",…

Change too early to say, but surgery was successful to sucessful

Show the code

df <- df %>% 
  mutate(outcome = fct_recode(outcome, "successful" = "too early to say, but surgery was successful"), 
         stage_of_develop = as.factor(stage_of_develop), 
         group = as.factor(group))

Donor sites

Show the code

table_donor <- df %>%
  separate(donor_tooth_recipient_region,
           into = c("donor", "recipient")) %>%
  mutate(donor2 = case_when(str_starts(donor, "1|2") ~ "maxilla",
                            TRUE ~ "mandibular")) %>%
  mutate(recipient2 = case_when(str_starts(recipient, "1|2") ~ "maxilla",
                                TRUE ~ "mandibular")) %>% 
  select(gender, group,donor, recipient ) %>% 
  janitor::tabyl(donor, recipient)

Show the code

table_donor %>% 
   gt()%>%
  tab_header(
    title = md("**Donor and recipient sites**")
  ) %>% 
   gt_theme_nytimes() %>% 
   tab_spanner(
    label = md("**Recipient site**"),
    columns = c("16":"47")
  ) %>% 
  gt_color_rows("16":"47")

Donor and recipient sites
donor	Recipient site
donor	16	26	35	36	37	46	47
18	2	2	0	2	1	2	0
28	2	4	1	7	1	1	0
38	1	0	0	7	4	1	0
48	0	0	0	3	1	3	1

Create a new var for donor sites

Show the code

df <- df %>%
  separate(donor_tooth_recipient_region,
           into = c("donor", "recipient")) %>%
  mutate(donor2 = case_when(str_starts(donor, "1|2") ~ "maxilla",
                            TRUE ~ "mandibular")) %>%
  mutate(recipient2 = case_when(str_starts(recipient, "1|2") ~ "maxilla",
                                TRUE ~ "mandibular")) %>%
  relocate(donor2, .before = donor) %>%
  relocate(recipient2, .after = donor2) %>% 
  unite("donor_recipient", donor2, recipient2, sep = "_") %>% 
  select(-c(donor, recipient, date_of_birth))

EDA

Overall

Show the code

df %>% 
  select(-patient_nr) %>% 
  describe_distribution() %>% 
  mutate(across(where(is.numeric), round, 1)) %>% 
  knitr::kable()

Variable	Mean	SD	IQR	Min	Max	Skewness	Kurtosis	n
age_surgery	17.1	1.8	2.0	13	22	0.4	0.6	46
total_surgery_time_min	54.9	18.8	21.2	25	120	1.4	3.2	46
extra_alveolar_time_s	59.9	46.3	30.0	10	180	1.6	1.8	46
donor_fitting_times	1.8	1.0	1.2	1	5	1.1	1.1	46

Age by surgery and gender

Show the code

df %>% 
  ggplot(aes(x = age_surgery, 
             fill = gender)) +
  geom_histogram(bins = 7) +
  facet_grid(group ~ gender)

Age by gender

Show the code

pp <- df %>% 
  ggplot(aes(y = age_surgery, 
             x = gender))  +
  geom_jitter(width = .2, alpha = .5) + 
  geom_boxplot(alpha = .5) +
  ylim(10, 25)

ggpval::add_pval(pp, pairs = list(c(1, 2)), test='t.test', alternative='two.sided')

Show the code

t.test(df$age_surgery ~ df$gender)


    Welch Two Sample t-test

data:  df$age_surgery by df$gender
t = -0.21, df = 42.451, p-value = 0.8347
alternative hypothesis: true difference in means between group Female and group Male is not equal to 0
95 percent confidence interval:
 -1.1402810  0.9252709
sample estimates:
mean in group Female   mean in group Male 
            17.06897             17.17647

Age by group

Show the code

pp <- df %>% 
  ggplot(aes(x = group, 
             y = age_surgery)) +
  geom_jitter(width = .2, alpha = .5) + 
  geom_boxplot(alpha = .5) +
  ylim(10, 25)

ggpval::add_pval(pp, pairs = list(c(1, 2)), test='t.test', alternative='two.sided')

Show the code

t.test(df$age_surgery ~ df$group)


    Welch Two Sample t-test

data:  df$age_surgery by df$group
t = 2.6531, df = 33.779, p-value = 0.01206
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
 0.3179813 2.4017157
sample estimates:
mean in group 0 mean in group 1 
       17.81818        16.45833

Do the groups have the same stage of development?

Show the code

mosaicplot(table(df$group, df$stage_of_develop), shade = T)

Whats is the outcome per group?

Show the code

mosaicplot(table(df$group, df$outcome), shade = T)

Surgery time per group

Show the code

pp <- df %>% 
  ggplot(aes(x = group, 
             y = total_surgery_time_min)) +
  geom_jitter(alpha = .5, width = .3) +
  geom_boxplot(alpha = .5) 

ggpval::add_pval(pp, pairs = list(c(1, 2)), test='t.test', alternative='two.sided')

Show the code

t.test(df$total_surgery_time_min ~ df$group)


    Welch Two Sample t-test

data:  df$total_surgery_time_min by df$group
t = 3.3737, df = 33.687, p-value = 0.001878
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
  6.834488 27.559452
sample estimates:
mean in group 0 mean in group 1 
       63.86364        46.66667

Extralveolar time per group

Show the code

pp <- df %>% 
  ggplot(aes(x = group, 
             y = extra_alveolar_time_s)) +
  geom_jitter(alpha = .5, width = .3) +
  geom_boxplot(alpha = .5) 

ggpval::add_pval(pp, pairs = list(c(1, 2)), test='t.test', alternative='two.sided')

Table 1

Show the code

df %>% 
  select(-patient_nr, 
         -date_of_surgery,
         -last_control) %>% 
  gtsummary::tbl_summary(by = group) %>% 
  gtsummary::bold_labels()

Characteristic	0, N = 22¹	1, N = 24¹
gender
Female	16 (73%)	13 (54%)
Male	6 (27%)	11 (46%)
age_surgery	17.50 (17.00, 19.00)	16.00 (16.00, 17.00)
total_surgery_time_min	60 (50, 70)	45 (40, 51)
extra_alveolar_time_s	55 (30, 82)	35 (30, 60)
donor_fitting_times
1	7 (32%)	15 (62%)
2	8 (36%)	5 (21%)
3	5 (23%)	4 (17%)
4	1 (4.5%)	0 (0%)
5	1 (4.5%)	0 (0%)
stage_of_develop
3	8 (36%)	9 (38%)
4	6 (27%)	12 (50%)
5	8 (36%)	3 (12%)
donor_recipient
mandibular_mandibular	10 (45%)	10 (42%)
mandibular_maxilla	1 (4.5%)	0 (0%)
maxilla_mandibular	8 (36%)	7 (29%)
maxilla_maxilla	3 (14%)	7 (29%)
outcome
extracted	1 (4.5%)	1 (4.2%)
successful	19 (86%)	22 (92%)
survived	2 (9.1%)	1 (4.2%)
¹ n (%); Median (IQR)

Table 2: Regression analysis

Is the time influenced by the group, adjusted by sex, age, surgery time, extraoral time, stage of development?

Model

Show the code

m1 <- glm(total_surgery_time_min ~ group + gender + extra_alveolar_time_s +
            donor_fitting_times + stage_of_develop + donor_recipient, 
          data = df)

Model evaluation

Show the code

performance::model_performance(m1)

# Indices of model performance

AIC     |     BIC |    R2 |   RMSE |  Sigma
-------------------------------------------
399.285 | 419.400 | 0.385 | 14.614 | 16.520

Show the code

check_heteroscedasticity(m1)

OK: Error variance appears to be homoscedastic (p = 0.332).

Regression table

Show the code

m1 %>% 
  gtsummary::tbl_regression() %>% 
  gtsummary::bold_labels() %>% 
  gtsummary::add_n(location = "level")

Characteristic	N	Beta	95% CI¹	p-value
group
0	22	—	—
1	24	-14	-25, -3.3	0.016
gender
Female	29	—	—
Male	17	3.7	-7.2, 15	0.5
extra_alveolar_time_s	46	-0.02	-0.18, 0.14	0.8
donor_fitting_times	46	3.5	-5.2, 12	0.4
stage_of_develop
3	17	—	—
4	18	3.5	-8.3, 15	0.6
5	11	4.1	-9.5, 18	0.6
donor_recipient
mandibular_mandibular	20	—	—
mandibular_maxilla	1	40	0.04, 79	0.058
maxilla_mandibular	15	-2.1	-14, 9.5	0.7
maxilla_maxilla	10	0.23	-13, 14	>0.9
¹ CI = Confidence Interval

Show the code

sjPlot::plot_model(m1, 
                   vline.color = "dark grey", 
                   sort.est = TRUE, 
                   title =  "Regression model for the effect of the 3D printing")

Show the code

report(m1)

We fitted a linear model (estimated using ML) to predict total_surgery_time_min
with group (formula: total_surgery_time_min ~ group + gender +
extra_alveolar_time_s + donor_fitting_times + stage_of_develop +
donor_recipient). The model's explanatory power is substantial (R2 = 0.39). The
model's intercept, corresponding to group = 0, is at 53.38 (95% CI [37.57,
69.20], t(36) = 6.61, p < .001). Within this model:

  - The effect of group [1] is statistically significant and negative (beta =
-14.36, 95% CI [-25.43, -3.28], t(36) = -2.54, p = 0.011; Std. beta = -0.76,
95% CI [-1.35, -0.17])
  - The effect of gender [Male] is statistically non-significant and positive
(beta = 3.73, 95% CI [-7.18, 14.65], t(36) = 0.67, p = 0.502; Std. beta = 0.20,
95% CI [-0.38, 0.78])
  - The effect of extra alveolar time s is statistically non-significant and
negative (beta = -0.02, 95% CI [-0.18, 0.14], t(36) = -0.27, p = 0.790; Std.
beta = -0.05, 95% CI [-0.45, 0.34])
  - The effect of donor fitting times is statistically non-significant and
positive (beta = 3.47, 95% CI [-5.25, 12.18], t(36) = 0.78, p = 0.435; Std.
beta = 0.18, 95% CI [-0.27, 0.63])
  - The effect of stage of develop [4] is statistically non-significant and
positive (beta = 3.52, 95% CI [-8.28, 15.32], t(36) = 0.59, p = 0.558; Std.
beta = 0.19, 95% CI [-0.44, 0.81])
  - The effect of stage of develop [5] is statistically non-significant and
positive (beta = 4.10, 95% CI [-9.47, 17.66], t(36) = 0.59, p = 0.554; Std.
beta = 0.22, 95% CI [-0.50, 0.94])
  - The effect of donor recipient [mandibular_maxilla] is statistically
significant and positive (beta = 39.71, 95% CI [0.04, 79.37], t(36) = 1.96, p =
0.050; Std. beta = 2.11, 95% CI [1.91e-03, 4.21])
  - The effect of donor recipient [maxilla_mandibular] is statistically
non-significant and negative (beta = -2.06, 95% CI [-13.58, 9.46], t(36) =
-0.35, p = 0.726; Std. beta = -0.11, 95% CI [-0.72, 0.50])
  - The effect of donor recipient [maxilla_maxilla] is statistically
non-significant and positive (beta = 0.23, 95% CI [-13.27, 13.73], t(36) =
0.03, p = 0.973; Std. beta = 0.01, 95% CI [-0.70, 0.73])

Standardized parameters were obtained by fitting the model on a standardized
version of the dataset. 95% Confidence Intervals (CIs) and p-values were
computed using a Wald t-distribution approximation., We fitted a linear model
(estimated using ML) to predict total_surgery_time_min with gender (formula:
total_surgery_time_min ~ group + gender + extra_alveolar_time_s +
donor_fitting_times + stage_of_develop + donor_recipient). The model's
explanatory power is substantial (R2 = 0.39). The model's intercept,
corresponding to gender = Female, is at 53.38 (95% CI [37.57, 69.20], t(36) =
6.61, p < .001). Within this model:

  - The effect of group [1] is statistically significant and negative (beta =
-14.36, 95% CI [-25.43, -3.28], t(36) = -2.54, p = 0.011; Std. beta = -0.76,
95% CI [-1.35, -0.17])
  - The effect of gender [Male] is statistically non-significant and positive
(beta = 3.73, 95% CI [-7.18, 14.65], t(36) = 0.67, p = 0.502; Std. beta = 0.20,
95% CI [-0.38, 0.78])
  - The effect of extra alveolar time s is statistically non-significant and
negative (beta = -0.02, 95% CI [-0.18, 0.14], t(36) = -0.27, p = 0.790; Std.
beta = -0.05, 95% CI [-0.45, 0.34])
  - The effect of donor fitting times is statistically non-significant and
positive (beta = 3.47, 95% CI [-5.25, 12.18], t(36) = 0.78, p = 0.435; Std.
beta = 0.18, 95% CI [-0.27, 0.63])
  - The effect of stage of develop [4] is statistically non-significant and
positive (beta = 3.52, 95% CI [-8.28, 15.32], t(36) = 0.59, p = 0.558; Std.
beta = 0.19, 95% CI [-0.44, 0.81])
  - The effect of stage of develop [5] is statistically non-significant and
positive (beta = 4.10, 95% CI [-9.47, 17.66], t(36) = 0.59, p = 0.554; Std.
beta = 0.22, 95% CI [-0.50, 0.94])
  - The effect of donor recipient [mandibular_maxilla] is statistically
significant and positive (beta = 39.71, 95% CI [0.04, 79.37], t(36) = 1.96, p =
0.050; Std. beta = 2.11, 95% CI [1.91e-03, 4.21])
  - The effect of donor recipient [maxilla_mandibular] is statistically
non-significant and negative (beta = -2.06, 95% CI [-13.58, 9.46], t(36) =
-0.35, p = 0.726; Std. beta = -0.11, 95% CI [-0.72, 0.50])
  - The effect of donor recipient [maxilla_maxilla] is statistically
non-significant and positive (beta = 0.23, 95% CI [-13.27, 13.73], t(36) =
0.03, p = 0.973; Std. beta = 0.01, 95% CI [-0.70, 0.73])

Standardized parameters were obtained by fitting the model on a standardized
version of the dataset. 95% Confidence Intervals (CIs) and p-values were
computed using a Wald t-distribution approximation., We fitted a linear model
(estimated using ML) to predict total_surgery_time_min with
extra_alveolar_time_s (formula: total_surgery_time_min ~ group + gender +
extra_alveolar_time_s + donor_fitting_times + stage_of_develop +
donor_recipient). The model's explanatory power is substantial (R2 = 0.39). The
model's intercept, corresponding to extra_alveolar_time_s = 0, is at 53.38 (95%
CI [37.57, 69.20], t(36) = 6.61, p < .001). Within this model:

  - The effect of group [1] is statistically significant and negative (beta =
-14.36, 95% CI [-25.43, -3.28], t(36) = -2.54, p = 0.011; Std. beta = -0.76,
95% CI [-1.35, -0.17])
  - The effect of gender [Male] is statistically non-significant and positive
(beta = 3.73, 95% CI [-7.18, 14.65], t(36) = 0.67, p = 0.502; Std. beta = 0.20,
95% CI [-0.38, 0.78])
  - The effect of extra alveolar time s is statistically non-significant and
negative (beta = -0.02, 95% CI [-0.18, 0.14], t(36) = -0.27, p = 0.790; Std.
beta = -0.05, 95% CI [-0.45, 0.34])
  - The effect of donor fitting times is statistically non-significant and
positive (beta = 3.47, 95% CI [-5.25, 12.18], t(36) = 0.78, p = 0.435; Std.
beta = 0.18, 95% CI [-0.27, 0.63])
  - The effect of stage of develop [4] is statistically non-significant and
positive (beta = 3.52, 95% CI [-8.28, 15.32], t(36) = 0.59, p = 0.558; Std.
beta = 0.19, 95% CI [-0.44, 0.81])
  - The effect of stage of develop [5] is statistically non-significant and
positive (beta = 4.10, 95% CI [-9.47, 17.66], t(36) = 0.59, p = 0.554; Std.
beta = 0.22, 95% CI [-0.50, 0.94])
  - The effect of donor recipient [mandibular_maxilla] is statistically
significant and positive (beta = 39.71, 95% CI [0.04, 79.37], t(36) = 1.96, p =
0.050; Std. beta = 2.11, 95% CI [1.91e-03, 4.21])
  - The effect of donor recipient [maxilla_mandibular] is statistically
non-significant and negative (beta = -2.06, 95% CI [-13.58, 9.46], t(36) =
-0.35, p = 0.726; Std. beta = -0.11, 95% CI [-0.72, 0.50])
  - The effect of donor recipient [maxilla_maxilla] is statistically
non-significant and positive (beta = 0.23, 95% CI [-13.27, 13.73], t(36) =
0.03, p = 0.973; Std. beta = 0.01, 95% CI [-0.70, 0.73])

Standardized parameters were obtained by fitting the model on a standardized
version of the dataset. 95% Confidence Intervals (CIs) and p-values were
computed using a Wald t-distribution approximation., We fitted a linear model
(estimated using ML) to predict total_surgery_time_min with donor_fitting_times
(formula: total_surgery_time_min ~ group + gender + extra_alveolar_time_s +
donor_fitting_times + stage_of_develop + donor_recipient). The model's
explanatory power is substantial (R2 = 0.39). The model's intercept,
corresponding to donor_fitting_times = 0, is at 53.38 (95% CI [37.57, 69.20],
t(36) = 6.61, p < .001). Within this model:

  - The effect of group [1] is statistically significant and negative (beta =
-14.36, 95% CI [-25.43, -3.28], t(36) = -2.54, p = 0.011; Std. beta = -0.76,
95% CI [-1.35, -0.17])
  - The effect of gender [Male] is statistically non-significant and positive
(beta = 3.73, 95% CI [-7.18, 14.65], t(36) = 0.67, p = 0.502; Std. beta = 0.20,
95% CI [-0.38, 0.78])
  - The effect of extra alveolar time s is statistically non-significant and
negative (beta = -0.02, 95% CI [-0.18, 0.14], t(36) = -0.27, p = 0.790; Std.
beta = -0.05, 95% CI [-0.45, 0.34])
  - The effect of donor fitting times is statistically non-significant and
positive (beta = 3.47, 95% CI [-5.25, 12.18], t(36) = 0.78, p = 0.435; Std.
beta = 0.18, 95% CI [-0.27, 0.63])
  - The effect of stage of develop [4] is statistically non-significant and
positive (beta = 3.52, 95% CI [-8.28, 15.32], t(36) = 0.59, p = 0.558; Std.
beta = 0.19, 95% CI [-0.44, 0.81])
  - The effect of stage of develop [5] is statistically non-significant and
positive (beta = 4.10, 95% CI [-9.47, 17.66], t(36) = 0.59, p = 0.554; Std.
beta = 0.22, 95% CI [-0.50, 0.94])
  - The effect of donor recipient [mandibular_maxilla] is statistically
significant and positive (beta = 39.71, 95% CI [0.04, 79.37], t(36) = 1.96, p =
0.050; Std. beta = 2.11, 95% CI [1.91e-03, 4.21])
  - The effect of donor recipient [maxilla_mandibular] is statistically
non-significant and negative (beta = -2.06, 95% CI [-13.58, 9.46], t(36) =
-0.35, p = 0.726; Std. beta = -0.11, 95% CI [-0.72, 0.50])
  - The effect of donor recipient [maxilla_maxilla] is statistically
non-significant and positive (beta = 0.23, 95% CI [-13.27, 13.73], t(36) =
0.03, p = 0.973; Std. beta = 0.01, 95% CI [-0.70, 0.73])

Standardized parameters were obtained by fitting the model on a standardized
version of the dataset. 95% Confidence Intervals (CIs) and p-values were
computed using a Wald t-distribution approximation., We fitted a linear model
(estimated using ML) to predict total_surgery_time_min with stage_of_develop
(formula: total_surgery_time_min ~ group + gender + extra_alveolar_time_s +
donor_fitting_times + stage_of_develop + donor_recipient). The model's
explanatory power is substantial (R2 = 0.39). The model's intercept,
corresponding to stage_of_develop = 3, is at 53.38 (95% CI [37.57, 69.20],
t(36) = 6.61, p < .001). Within this model:

  - The effect of group [1] is statistically significant and negative (beta =
-14.36, 95% CI [-25.43, -3.28], t(36) = -2.54, p = 0.011; Std. beta = -0.76,
95% CI [-1.35, -0.17])
  - The effect of gender [Male] is statistically non-significant and positive
(beta = 3.73, 95% CI [-7.18, 14.65], t(36) = 0.67, p = 0.502; Std. beta = 0.20,
95% CI [-0.38, 0.78])
  - The effect of extra alveolar time s is statistically non-significant and
negative (beta = -0.02, 95% CI [-0.18, 0.14], t(36) = -0.27, p = 0.790; Std.
beta = -0.05, 95% CI [-0.45, 0.34])
  - The effect of donor fitting times is statistically non-significant and
positive (beta = 3.47, 95% CI [-5.25, 12.18], t(36) = 0.78, p = 0.435; Std.
beta = 0.18, 95% CI [-0.27, 0.63])
  - The effect of stage of develop [4] is statistically non-significant and
positive (beta = 3.52, 95% CI [-8.28, 15.32], t(36) = 0.59, p = 0.558; Std.
beta = 0.19, 95% CI [-0.44, 0.81])
  - The effect of stage of develop [5] is statistically non-significant and
positive (beta = 4.10, 95% CI [-9.47, 17.66], t(36) = 0.59, p = 0.554; Std.
beta = 0.22, 95% CI [-0.50, 0.94])
  - The effect of donor recipient [mandibular_maxilla] is statistically
significant and positive (beta = 39.71, 95% CI [0.04, 79.37], t(36) = 1.96, p =
0.050; Std. beta = 2.11, 95% CI [1.91e-03, 4.21])
  - The effect of donor recipient [maxilla_mandibular] is statistically
non-significant and negative (beta = -2.06, 95% CI [-13.58, 9.46], t(36) =
-0.35, p = 0.726; Std. beta = -0.11, 95% CI [-0.72, 0.50])
  - The effect of donor recipient [maxilla_maxilla] is statistically
non-significant and positive (beta = 0.23, 95% CI [-13.27, 13.73], t(36) =
0.03, p = 0.973; Std. beta = 0.01, 95% CI [-0.70, 0.73])

Standardized parameters were obtained by fitting the model on a standardized
version of the dataset. 95% Confidence Intervals (CIs) and p-values were
computed using a Wald t-distribution approximation. and We fitted a linear
model (estimated using ML) to predict total_surgery_time_min with
donor_recipient (formula: total_surgery_time_min ~ group + gender +
extra_alveolar_time_s + donor_fitting_times + stage_of_develop +
donor_recipient). The model's explanatory power is substantial (R2 = 0.39). The
model's intercept, corresponding to donor_recipient = mandibular_mandibular, is
at 53.38 (95% CI [37.57, 69.20], t(36) = 6.61, p < .001). Within this model:

  - The effect of group [1] is statistically significant and negative (beta =
-14.36, 95% CI [-25.43, -3.28], t(36) = -2.54, p = 0.011; Std. beta = -0.76,
95% CI [-1.35, -0.17])
  - The effect of gender [Male] is statistically non-significant and positive
(beta = 3.73, 95% CI [-7.18, 14.65], t(36) = 0.67, p = 0.502; Std. beta = 0.20,
95% CI [-0.38, 0.78])
  - The effect of extra alveolar time s is statistically non-significant and
negative (beta = -0.02, 95% CI [-0.18, 0.14], t(36) = -0.27, p = 0.790; Std.
beta = -0.05, 95% CI [-0.45, 0.34])
  - The effect of donor fitting times is statistically non-significant and
positive (beta = 3.47, 95% CI [-5.25, 12.18], t(36) = 0.78, p = 0.435; Std.
beta = 0.18, 95% CI [-0.27, 0.63])
  - The effect of stage of develop [4] is statistically non-significant and
positive (beta = 3.52, 95% CI [-8.28, 15.32], t(36) = 0.59, p = 0.558; Std.
beta = 0.19, 95% CI [-0.44, 0.81])
  - The effect of stage of develop [5] is statistically non-significant and
positive (beta = 4.10, 95% CI [-9.47, 17.66], t(36) = 0.59, p = 0.554; Std.
beta = 0.22, 95% CI [-0.50, 0.94])
  - The effect of donor recipient [mandibular_maxilla] is statistically
significant and positive (beta = 39.71, 95% CI [0.04, 79.37], t(36) = 1.96, p =
0.050; Std. beta = 2.11, 95% CI [1.91e-03, 4.21])
  - The effect of donor recipient [maxilla_mandibular] is statistically
non-significant and negative (beta = -2.06, 95% CI [-13.58, 9.46], t(36) =
-0.35, p = 0.726; Std. beta = -0.11, 95% CI [-0.72, 0.50])
  - The effect of donor recipient [maxilla_maxilla] is statistically
non-significant and positive (beta = 0.23, 95% CI [-13.27, 13.73], t(36) =
0.03, p = 0.973; Std. beta = 0.01, 95% CI [-0.70, 0.73])

Standardized parameters were obtained by fitting the model on a standardized
version of the dataset. 95% Confidence Intervals (CIs) and p-values were
computed using a Wald t-distribution approximation.

Show the code

check_normality(m1) %>% plot(type = "qq")

Show the code

check_heteroscedasticity(m1) %>% plot()

Nonparametric bootstrap regression

1. Bootstrap the data with resampling

Show the code

set.seed(123)
boot_data <- bootstraps(df, times = 1000)

2. Produce 1000 models

Show the code

boot_models <- boot_data %>%
  mutate(model = map(
    splits,
    ~ glm(
      total_surgery_time_min ~ group + gender + extra_alveolar_time_s +
        donor_fitting_times + stage_of_develop + donor_recipient,
      data = .
    )
  ),
  coefs = map(model, tidy, conf.int = TRUE)) # optional "conf.int = TRUE"

3. Plot estimates distribution of 1000 estimates

Show the code

boot_coefs <- boot_models %>%
  unnest(coefs) %>% 
  select(term, estimate)

Show the code

boot_coefs %>%
  ggplot(aes(x = estimate)) +
  geom_histogram() +
  facet_wrap(~term, scales = "free")

4. Get and visualize predictions

Show the code

# get predictions
boot_aug <- boot_models %>% 
  mutate(augmented = map(model, augment)) %>% 
  unnest(augmented) %>% 
  select(-splits, -model)

# get median estimates, median 95% CIs
nonpar_med_boot_preds <- boot_aug %>% 
  group_by(group) %>% 
  summarise(
    predicted = median(.fitted ),
    conf.low  = quantile(.fitted, 0.025),
    conf.high = quantile(.fitted, 0.975)) %>% 
  mutate(model = "bootstrapped") %>% 
  select(model, everything())

nonpar_med_boot_preds

# A tibble: 2 × 5
  model        group predicted conf.low conf.high
  <chr>        <fct>     <dbl>    <dbl>     <dbl>
1 bootstrapped 0          62.1     46.9     110  
2 bootstrapped 1          46.4     33.8      61.9

Compare with the traditional model

Show the code

par_avg_preds <- ggeffects::ggeffect(m1, terms = "group") %>% 
  tibble() %>% 
  mutate(model = "linear") %>% 
  select(model, group = x, predicted, conf.low, conf.high)
par_avg_preds

# A tibble: 2 × 5
  model  group predicted conf.low conf.high
  <chr>  <fct>     <dbl>    <dbl>     <dbl>
1 linear 0          62.4     54.6      70.1
2 linear 1          48.0     40.6      55.4

Show the code

preds <- rbind(nonpar_med_boot_preds, par_avg_preds)

Show the code

# plot row data and predicted CIs from both models
ggplot() +
  geom_jitter(data = df, aes(group, total_surgery_time_min), 
              width = .2, alpha = 0.2, size = 2)+
  geom_pointrange(data = preds, aes(x = group, y = predicted, 
                  ymin = conf.low, ymax = conf.high, color = model),
                  position=position_dodge(width=.6), size = 1)

Check with a robust model

Show the code

rm <- lmrob(total_surgery_time_min ~ group + gender + extra_alveolar_time_s +
            donor_fitting_times + stage_of_develop + donor_recipient, 
          data = df)

rob_avg_preds <- ggeffects::ggeffect(rm, terms = "group") %>% 
  tibble() %>% 
  mutate(model = "robust") %>% 
  select(model, group = x, predicted, conf.low, conf.high)

# combine two datasets
preds <- rbind(nonpar_med_boot_preds, par_avg_preds, rob_avg_preds)

ggplot() +
  geom_jitter(
    data = df,
    aes(group, total_surgery_time_min),
    width = .2,
    alpha = 0.2,
    size = 2
  ) +
  geom_pointrange(
    data = preds,
    aes(
      x = group,
      y = predicted,
      ymin = conf.low,
      ymax = conf.high,
      color = model
    ),
    position = position_dodge(width = .6),
    size = 1
  )

Summary

Conventional regression (glm) shows the effect of 3D use over time. The bootstrap model minimizes the effect, due to the small sample size.

If we consider that surgical time is a variable that has an impact on the survival of autotransplanted teeth, then the use of 3D models could contribute to improve the results of this biological and conservative intervention.