Burnout

library(tidyverse)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.0     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.0
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(gt)
library(gtsummary)

Demographics

load(file="data_clean.RData")
base<-base %>%filter(origin=="Fleni")
base %>% select(age, gender, family_income, personal_income,only_provider, laboral, work_type, work_md) %>%
  tbl_summary(by=work_type) %>% 
  add_overall %>%
  add_p() %>%
  add_q() %>%
  modify_caption("**Table 1** Demographics features for center (N = {N})")
Table 1 Demographics features for center (N = 615)
Characteristic Overall, N = 6151 Public-facing customer service (receptionists and secretaries), N = 781 Direct patient care (physician, nurse, physiotherapist, therapist, etc.), N = 3691 Indirect patient care (biochemist, imaging studies report, etc.), N = 711 ** Administrative tasks related to healthcare without public-facing duties (commercial administrative sector)**, N = 971 p-value2 q-value3
age 40 (33, 51) 37 (31, 46) 42 (34, 52) 41 (33, 50) 40 (32, 46) 0.032 0.032
    Unknown 348 40 220 32 56

gender




0.005 0.008
    Male 190 (31%) 14 (18%) 107 (29%) 31 (44%) 38 (39%)

    Female 423 (69%) 64 (82%) 260 (70%) 40 (56%) 59 (61%)

    Other 2 (0.3%) 0 (0%) 2 (0.5%) 0 (0%) 0 (0%)

family_income






    Less than $100.000 5 (0.8%) 1 (1.3%) 3 (0.8%) 0 (0%) 1 (1.0%)

    $100.001 to $ 250.000 104 (17%) 29 (37%) 43 (12%) 7 (9.9%) 25 (26%)

    $250.001 to 400.000 177 (29%) 25 (32%) 99 (27%) 22 (31%) 31 (32%)

    More than $400.000 329 (53%) 23 (29%) 224 (61%) 42 (59%) 40 (41%)

personal_income






    Less than $100.000 9 (1.5%) 3 (3.8%) 3 (0.8%) 0 (0%) 3 (3.1%)

    $100.001 to $ 250.000 196 (32%) 62 (79%) 75 (20%) 15 (21%) 44 (45%)

    $250.001 to 400.000 225 (37%) 12 (15%) 157 (43%) 30 (42%) 26 (27%)

    More than $400.000 185 (30%) 1 (1.3%) 134 (36%) 26 (37%) 24 (25%)

only_provider 396 (64%) 37 (47%) 251 (68%) 48 (68%) 60 (62%) 0.006 0.008
laboral






    Private practice 519 (85%) 67 (86%) 303 (82%) 57 (81%) 92 (95%)

    Public 15 (2.4%) 6 (7.7%) 6 (1.6%) 2 (2.9%) 1 (1.0%)

    Both mostly public 18 (2.9%) 3 (3.8%) 12 (3.3%) 3 (4.3%) 0 (0%)

    Both mostly private 62 (10%) 2 (2.6%) 48 (13%) 8 (11%) 4 (4.1%)

    Unknown 1 0 0 1 0

work_md 189 (31%) 0 (0%) 162 (44%) 24 (34%) 3 (3.1%) <0.001 <0.001
1 Median (IQR); n (%)
2 Kruskal-Wallis rank sum test; Fisher’s exact test; Pearson’s Chi-squared test
3 False discovery rate correction for multiple testing 
base %>% select(age, gender, family_income, personal_income,only_provider, laboral, work_type, work_md) %>%
  tbl_summary(by=work_md) %>% 
  add_p() %>%
  add_q() %>%
  modify_caption("**Table 2** Demographics MD vs non-MD (N = {N})")
Table 2 Demographics MD vs non-MD (N = 615)
Characteristic No, N = 4261 Yes, N = 1891 p-value2 q-value3
age 41 (33, 51) 39 (34, 50) >0.9 >0.9
    Unknown 244 104

gender

0.002 0.003
    Male 115 (27%) 75 (40%)

    Female 310 (73%) 113 (60%)

    Other 1 (0.2%) 1 (0.5%)

family_income

<0.001 <0.001
    Less than $100.000 5 (1.2%) 0 (0%)

    $100.001 to $ 250.000 98 (23%) 6 (3.2%)

    $250.001 to 400.000 131 (31%) 46 (24%)

    More than $400.000 192 (45%) 137 (72%)

personal_income

<0.001 <0.001
    Less than $100.000 8 (1.9%) 1 (0.5%)

    $100.001 to $ 250.000 182 (43%) 14 (7.4%)

    $250.001 to 400.000 157 (37%) 68 (36%)

    More than $400.000 79 (19%) 106 (56%)

only_provider 248 (58%) 148 (78%) <0.001 <0.001
laboral

0.010 0.012
    Private practice 367 (86%) 152 (81%)

    Public 14 (3.3%) 1 (0.5%)

    Both mostly public 10 (2.3%) 8 (4.3%)

    Both mostly private 35 (8.2%) 27 (14%)

    Unknown 0 1

work_type

<0.001 <0.001
    Public-facing customer service (receptionists and secretaries) 78 (18%) 0 (0%)

    Direct patient care (physician, nurse, physiotherapist, therapist, etc.) 207 (49%) 162 (86%)

    Indirect patient care (biochemist, imaging studies report, etc.) 47 (11%) 24 (13%)

     Administrative tasks related to healthcare without public-facing duties (commercial administrative sector) 94 (22%) 3 (1.6%)

1 Median (IQR); n (%)
2 Wilcoxon rank sum test; Fisher’s exact test; Pearson’s Chi-squared test
3 False discovery rate correction for multiple testing 
base %>% 
  filter(work_md=="Yes") %>%
  select(age, 
         gender, years_of_graduation,
         work_type,
         family_income, 
         personal_income,only_provider, 
         laboral,   
         laboral_dichot, work_hours,
         work_number,
         work_teaching,
         work_emergency,
         work_uti) %>%
  tbl_summary(by=laboral_dichot) %>% 
  add_overall()%>%
  add_p() %>%
  add_q() %>%
  modify_caption("**Table 3** Demographics of MD (N = {N})")
Table 3 Demographics of MD (N = 188)
Characteristic Overall, N = 1881 Private, N = 1791 Public, N = 91 p-value2 q-value3
age 39 (34, 50) 40 (33, 51) 37 (37, 42) >0.9 >0.9
    Unknown 103 97 6

gender


0.2 0.7
    Male 75 (40%) 69 (39%) 6 (67%)

    Female 112 (60%) 109 (61%) 3 (33%)

    Other 1 (0.5%) 1 (0.6%) 0 (0%)

years_of_graduation 13 (7, 23) 13 (7, 23) 12 (6, 31) >0.9 >0.9
work_type


0.7 >0.9
    Public-facing customer service (receptionists and secretaries) 0 (0%) 0 (0%) 0 (0%)

    Direct patient care (physician, nurse, physiotherapist, therapist, etc.) 162 (86%) 153 (85%) 9 (100%)

    Indirect patient care (biochemist, imaging studies report, etc.) 23 (12%) 23 (13%) 0 (0%)

     Administrative tasks related to healthcare without public-facing duties (commercial administrative sector) 3 (1.6%) 3 (1.7%) 0 (0%)

family_income


>0.9 >0.9
    Less than $100.000 0 (0%) 0 (0%) 0 (0%)

    $100.001 to $ 250.000 6 (3.2%) 6 (3.4%) 0 (0%)

    $250.001 to 400.000 45 (24%) 43 (24%) 2 (22%)

    More than $400.000 137 (73%) 130 (73%) 7 (78%)

personal_income


>0.9 >0.9
    Less than $100.000 1 (0.5%) 1 (0.6%) 0 (0%)

    $100.001 to $ 250.000 14 (7.4%) 14 (7.8%) 0 (0%)

    $250.001 to 400.000 67 (36%) 64 (36%) 3 (33%)

    More than $400.000 106 (56%) 100 (56%) 6 (67%)

only_provider 147 (78%) 139 (78%) 8 (89%) 0.7 >0.9
laboral


<0.001 <0.001
    Private practice 152 (81%) 152 (85%) 0 (0%)

    Public 1 (0.5%) 0 (0%) 1 (11%)

    Both mostly public 8 (4.3%) 0 (0%) 8 (89%)

    Both mostly private 27 (14%) 27 (15%) 0 (0%)

work_hours 48 (40, 60) 47 (40, 59) 50 (40, 60) >0.9 >0.9
work_number


0.2 0.7
    0 2 (1.1%) 2 (1.1%) 0 (0%)

    1 51 (27%) 51 (28%) 0 (0%)

    2 52 (28%) 50 (28%) 2 (22%)

    3 60 (32%) 55 (31%) 5 (56%)

    4 15 (8.0%) 13 (7.3%) 2 (22%)

    5 6 (3.2%) 6 (3.4%) 0 (0%)

    6 2 (1.1%) 2 (1.1%) 0 (0%)

work_teaching 71 (38%) 65 (36%) 6 (67%) 0.084 0.5
work_emergency 59 (31%) 57 (32%) 2 (22%) 0.7 >0.9
work_uti 22 (12%) 21 (12%) 1 (11%) >0.9 >0.9
1 Median (IQR); n (%)
2 Wilcoxon rank sum test; Fisher’s exact test
3 False discovery rate correction for multiple testing

Neuro-Psychiatric scores

base %>% select(dass_stress, dass_stress_dx,dass_anxiety,dass_anxiety_dx,  dass_depression, dass_depression_dx, work_md)%>%
  tbl_summary(by=work_md)%>% 
  add_overall()%>%
  add_p() %>%
  add_q() %>%
  modify_caption("**Table 4** Depresion/Anxiety/Strees in MD vs non-MD (N = {N})")
Table 4 Depresion/Anxiety/Strees in MD vs non-MD (N = 615)
Characteristic Overall, N = 6151 No, N = 4261 Yes, N = 1891 p-value2 q-value3
dass_stress 12 (6, 20) 12 (6, 18) 14 (8, 22) 0.010 0.062
    Unknown 250 152 98

dass_stress_dx


0.2 0.6
    Normal 116 (35%) 91 (38%) 25 (29%)

    Mild 119 (36%) 90 (37%) 29 (34%)

    Moderate 71 (22%) 48 (20%) 23 (27%)

    Severe 16 (4.9%) 10 (4.1%) 6 (7.0%)

    Extremely severe 6 (1.8%) 3 (1.2%) 3 (3.5%)

    Unknown 287 184 103

dass_anxiety 4 (0, 8) 4 (0, 8) 4 (0, 8) 0.5 0.8
    Unknown 147 92 55

dass_anxiety_dx


0.7 0.8
    Normal 178 (55%) 125 (54%) 53 (59%)

    Mild 32 (10.0%) 25 (11%) 7 (7.8%)

    Moderate 52 (16%) 38 (16%) 14 (16%)

    Severe 23 (7.2%) 15 (6.5%) 8 (8.9%)

    Extremely severe 36 (11%) 28 (12%) 8 (8.9%)

    Unknown 294 195 99

dass_depression 6 (0, 10) 6 (0, 10) 6 (2, 12) 0.7 0.8
    Unknown 177 106 71

dass_depression_dx


0.8 0.8
    Normal 187 (57%) 137 (58%) 50 (56%)

    Mild 50 (15%) 34 (14%) 16 (18%)

    Moderate 64 (20%) 48 (20%) 16 (18%)

    Severe 17 (5.2%) 11 (4.7%) 6 (6.7%)

    Extremely severe 8 (2.5%) 6 (2.5%) 2 (2.2%)

    Unknown 289 190 99

1 Median (IQR); n (%)
2 Wilcoxon rank sum test; Fisher’s exact test; Pearson’s Chi-squared test
3 False discovery rate correction for multiple testing 
library("plot3D")
Warning: package 'plot3D' was built under R version 4.3.3
base$md<-as.integer(base$work_md=="Yes")
scatter3D(x=base$dass_stress, y=base$dass_anxiety, z=base$dass_depression,
          pch = 18, alpha=.5, cex = 0.7,phi = 20,
          col.var = base$md, 
          col = c("#1B9E77", "#D95F02"),
          ticktype = "detailed",
          colkey = list(at = c(3,35), side = 1,
          addlines = TRUE, length = 0.2, width = 0.2,
          labels = c("non-MD", "MD")),
          main = "DASS MD vs nonMD", xlab = "Stress",
          ylab ="Anxiety", zlab = "Depression")

base%>%
  ggplot(aes(x=dass_anxiety, y=dass_depression, color=work_md))+
  geom_jitter(width = .5, height = .5, alpha=.5)+scale_fill_manual(values=c("#1B9E77", "#D95F02"))+
  geom_hline(yintercept=9, linetype="dashed")+
  geom_hline(yintercept=12, linetype="dashed")+
  geom_hline(yintercept=20, linetype="dashed")+
  geom_hline(yintercept=27, linetype="dashed")+
  geom_vline(xintercept=6, linetype="dashed")+
  geom_vline(xintercept=9, linetype="dashed")+
  geom_vline(xintercept=14, linetype="dashed")+
  geom_vline(xintercept=19, linetype="dashed")+
  theme_minimal()+
  ggtitle("DASS MD vs nonMD")

Maslach

base %>% select(maslach_emotional_exhaustion, maslach_depersonalization,maslach_personal_accomplishment, work_md)%>%
  tbl_summary(by=work_md)%>% 
  add_overall()%>%
  add_p() %>%
  add_q() %>%
  modify_caption("**Table 5** Raw Maslach dimensions in MD vs non-MD (N = {N})")
Table 5 Raw Maslach dimensions in MD vs non-MD (N = 615)
Characteristic Overall, N = 6151 No, N = 4261 Yes, N = 1891 p-value2 q-value3
maslach_emotional_exhaustion 25 (16, 36) 22 (15, 32) 33 (22, 40) <0.001 <0.001
    Unknown 2 2 0

maslach_depersonalization 7.0 (5.0, 11.0) 7.0 (5.0, 10.0) 9.0 (6.0, 12.0) <0.001 <0.001
    Unknown 1 1 0

maslach_personal_accomplishment 29 (24, 34) 30 (24, 34) 29 (24, 33) 0.4 0.4
1 Median (IQR)
2 Wilcoxon rank sum test
3 False discovery rate correction for multiple testing 
base %>% select(maslach_emotional_exhaustion_index, maslach_depersonalization_index,maslach_personal_accomplishment_index, work_md)%>%
  tbl_summary(by=work_md)%>% 
  add_overall()%>%
  add_p() %>%
  add_q() %>%
  modify_caption("**Table 6** adjusted Maslach dimensions in MD vs non-MD (N = {N})")
Table 6 adjusted Maslach dimensions in MD vs non-MD (N = 615)
Characteristic Overall, N = 6151 No, N = 4261 Yes, N = 1891 p-value2 q-value3
maslach_emotional_exhaustion_index 0.46 (0.30, 0.67) 0.41 (0.28, 0.59) 0.61 (0.41, 0.74) <0.001 <0.001
    Unknown 2 2 0

maslach_depersonalization_index 0.23 (0.17, 0.37) 0.23 (0.17, 0.33) 0.30 (0.20, 0.40) <0.001 <0.001
    Unknown 1 1 0

maslach_personal_accomplishment_index 0.69 (0.57, 0.81) 0.71 (0.57, 0.81) 0.69 (0.57, 0.79) 0.4 0.4
1 Median (IQR)
2 Wilcoxon rank sum test
3 False discovery rate correction for multiple testing

Correlation Neuropsychiatric vs Maslach

library(corrplot)
corrplot 0.92 loaded
a<-base %>% select(maslach_personal_accomplishment_index,
                maslach_depersonalization_index,
                maslach_emotional_exhaustion_index,
                dass_anxiety_index,
                dass_depression_index,
                dass_stress_index)

a<-scale(a)
a<-na.omit(a)
M<-cor(a)

cor.mtest <- function(mat, ...) {
    mat <- as.matrix(mat)
    n <- ncol(mat)
    p.mat<- matrix(NA, n, n)
    diag(p.mat) <- 0
    for (i in 1:(n - 1)) {
        for (j in (i + 1):n) {
            tmp <- cor.test(mat[, i], mat[, j], ...)
            p.mat[i, j] <- p.mat[j, i] <- tmp$p.value
        }
    }
  colnames(p.mat) <- rownames(p.mat) <- colnames(mat)
  p.mat
}
# matrix of the p-value of the correlation
p.mat <- cor.mtest(a)

corrplot(M, method="color",  
         type="upper", order="hclust", 
         addCoef.col = "black", # Add coefficient of correlation
         tl.col="black", #Text label color and rotation
         # Combine with significance
         p.mat = p.mat, sig.level = 0.01, 
         # hide correlation coefficient on the principal diagonal
         diag=FALSE)

Modelos explicativos del stress

library(psych)

Attaching package: 'psych'
The following objects are masked from 'package:ggplot2':

    %+%, alpha

Emotional exhaustion explanatory model

Model 1

base_model<-base %>% select(maslach_emotional_exhaustion_index,       
                            maslach_personal_accomplishment_index,
                            maslach_depersonalization_index, work_md,
                            work_hours, personal_income, only_provider,
                            work_number) %>%   
  rename(emotional_e=maslach_emotional_exhaustion_index,
         personal_a=maslach_personal_accomplishment_index,
         depers=maslach_depersonalization_index,
         md=work_md)  

base_model$md<-as.numeric(base_model$md) 
base_model$only_provider<-as.numeric(base_model$only_provider)    

mediate(emotional_e+personal_a+depers~md+work_hours+only_provider, data = base_model, n.iter = 10000) %>%   print(large)


Mediation/Moderation Analysis 
Call: mediate(y = emotional_e + personal_a + depers ~ md + work_hours + 
    only_provider, data = base_model, n.iter = 10000)

The DV (Y) was  emotional_e personal_a depers . The IV (X) was  md work_hours only_provider . The mediating variable(s) =  .Call: mediate(y = emotional_e + personal_a + depers ~ md + work_hours + 
    only_provider, data = base_model, n.iter = 10000)

No mediator specified leads to traditional regression 
              emotional_e   se     t  df     Prob
Intercept            0.25 0.03  7.52 611 1.93e-13
md                   0.13 0.02  6.27 611 6.70e-10
work_hours           0.00 0.00  3.79 611 1.66e-04
only_provider       -0.04 0.02 -2.20 611 2.83e-02

R = 0.31 R2 = 0.1   F = 22.02 on 3 and 611 DF   p-value:  1.49e-13 
              personal_a   se     t  df      Prob
Intercept           0.69 0.02 27.88 611 5.84e-111
md                 -0.01 0.02 -0.83 611  4.07e-01
work_hours          0.00 0.00 -0.87 611  3.87e-01
only_provider       0.04 0.01  2.93 611  3.56e-03

R = 0.12 R2 = 0.01   F = 3.04 on 3 and 611 DF   p-value:  0.0285 
              depers   se    t  df     Prob
Intercept       0.19 0.02 8.12 611 2.50e-15
md              0.04 0.01 3.08 611 2.15e-03
work_hours      0.00 0.00 2.06 611 4.00e-02
only_provider   0.00 0.01 0.33 611 7.41e-01

R = 0.17 R2 = 0.03   F = 6.18 on 3 and 611 DF   p-value:  0.000387

Model 2

base_model<-base %>% select(maslach_emotional_exhaustion_index,
                            maslach_personal_accomplishment_index, maslach_depersonalization_index,
                      dass_depression_index,
                      dass_anxiety_index, dass_stress_index, work_md, work_hours, personal_income, only_provider, work_number) %>%
  rename(emotional_e=maslach_emotional_exhaustion_index,
         personal_a=maslach_personal_accomplishment_index,
         depers=maslach_depersonalization_index,
         depress=dass_depression_index,
         anx=dass_anxiety_index, stress=dass_stress_index, 
         md=work_md)

base_model$md<-as.numeric(base_model$md)
base_model$only_provider<-as.numeric(base_model$only_provider)



mediate(emotional_e+personal_a+depers~ 
          depress+
          anx +stress+md+
          work_hours+
          only_provider, data = base_model, n.iter = 10000) %>%
  print(large)


Mediation/Moderation Analysis 
Call: mediate(y = emotional_e + personal_a + depers ~ depress + anx + 
    stress + md + work_hours + only_provider, data = base_model, 
    n.iter = 10000)

The DV (Y) was  emotional_e personal_a depers . The IV (X) was  depress anx stress md work_hours only_provider . The mediating variable(s) =  .Call: mediate(y = emotional_e + personal_a + depers ~ depress + anx + 
    stress + md + work_hours + only_provider, data = base_model, 
    n.iter = 10000)

No mediator specified leads to traditional regression 
              emotional_e   se     t  df     Prob
Intercept            0.12 0.03  4.13 608 4.11e-05
depress              0.32 0.05  5.90 608 5.89e-09
anx                  0.13 0.05  2.33 608 1.99e-02
stress               0.31 0.05  6.27 608 6.89e-10
md                   0.11 0.02  6.34 608 4.44e-10
work_hours           0.00 0.00  3.07 608 2.22e-03
only_provider       -0.02 0.02 -1.17 608 2.43e-01

R = 0.61 R2 = 0.37   F = 59.54 on 6 and 608 DF   p-value:  6.18e-58 
              personal_a   se     t  df      Prob
Intercept           0.73 0.02 29.40 608 7.59e-119
depress            -0.27 0.05 -5.93 608  5.12e-09
anx                 0.02 0.05  0.36 608  7.22e-01
stress              0.02 0.04  0.39 608  6.94e-01
md                 -0.01 0.01 -0.68 608  4.95e-01
work_hours          0.00 0.00 -0.62 608  5.38e-01
only_provider       0.03 0.01  2.12 608  3.43e-02

R = 0.3 R2 = 0.09   F = 10.19 on 6 and 608 DF   p-value:  9.28e-11 
              depers   se    t  df     Prob
Intercept       0.15 0.02 6.48 608 1.86e-10
depress         0.10 0.04 2.19 608 2.88e-02
anx             0.03 0.04 0.76 608 4.49e-01
stress          0.07 0.04 1.84 608 6.65e-02
md              0.04 0.01 2.78 608 5.55e-03
work_hours      0.00 0.00 1.64 608 1.02e-01
only_provider   0.01 0.01 0.83 608 4.08e-01

R = 0.27 R2 = 0.07   F = 7.98 on 6 and 608 DF   p-value:  2.65e-08 
#  model <- psych::lmCor(TC + TS ~ BC + BS, data = dancer)
# > summary(model)

Model 3

base_model<-base %>% select(maslach_emotional_exhaustion_index,
                            maslach_personal_accomplishment_index, maslach_depersonalization_index,
                      dass_depression_index,
                      dass_anxiety_index, dass_stress_index, work_md, work_hours, personal_income, only_provider, work_number) %>%
  rename(emotional_e=maslach_emotional_exhaustion_index,
         personal_a=maslach_personal_accomplishment_index,
         depers=maslach_depersonalization_index,
         depress=dass_depression_index,
         anx=dass_anxiety_index, stress=dass_stress_index, 
         md=work_md)

base_model$md<-as.numeric(base_model$md)
base_model$only_provider<-as.numeric(base_model$only_provider)

base_model$personal_income<-as.numeric(base_model$personal_income)


mediate(emotional_e+depers~ 
          depress+
          anx +stress+md+
          work_hours+
          only_provider+ personal_income, data = base_model, n.iter = 10000) %>%
  print(large)


Mediation/Moderation Analysis 
Call: mediate(y = emotional_e + depers ~ depress + anx + stress + md + 
    work_hours + only_provider + personal_income, data = base_model, 
    n.iter = 10000)

The DV (Y) was  emotional_e depers . The IV (X) was  depress anx stress md work_hours only_provider personal_income . The mediating variable(s) =  .Call: mediate(y = emotional_e + depers ~ depress + anx + stress + md + 
    work_hours + only_provider + personal_income, data = base_model, 
    n.iter = 10000)

No mediator specified leads to traditional regression 
                emotional_e   se     t  df     Prob
Intercept              0.16 0.04  4.46 607 9.99e-06
depress                0.32 0.05  5.84 607 8.56e-09
anx                    0.11 0.05  2.01 607 4.51e-02
stress                 0.31 0.05  6.37 607 3.66e-10
md                     0.12 0.02  6.59 607 9.56e-11
work_hours             0.00 0.00  3.06 607 2.32e-03
only_provider         -0.01 0.02 -0.43 607 6.65e-01
personal_income       -0.02 0.01 -1.86 607 6.31e-02

R = 0.61 R2 = 0.37   F = 51.73 on 7 and 607 DF   p-value:  9.07e-58 
                depers   se     t  df     Prob
Intercept         0.16 0.03  5.42 607 8.78e-08
depress           0.10 0.04  2.18 607 2.96e-02
anx               0.03 0.04  0.71 607 4.77e-01
stress            0.07 0.04  1.85 607 6.53e-02
md                0.04 0.02  2.66 607 7.98e-03
work_hours        0.00 0.00  1.64 607 1.03e-01
only_provider     0.01 0.01  0.85 607 3.95e-01
personal_income   0.00 0.01 -0.22 607 8.26e-01

R = 0.27 R2 = 0.07   F = 6.84 on 7 and 607 DF   p-value:  7.8e-08

Model 4

base_model<-base %>% select(maslach_emotional_exhaustion_index,
                            maslach_personal_accomplishment_index, maslach_depersonalization_index,
                      dass_depression_index,
                      dass_anxiety_index, dass_stress_index, work_md, work_hours, personal_income, only_provider, work_number) %>%
  rename(emotional_e=maslach_emotional_exhaustion_index,
         personal_a=maslach_personal_accomplishment_index,
         depers=maslach_depersonalization_index,
         depress=dass_depression_index,
         anx=dass_anxiety_index, stress=dass_stress_index, 
         md=work_md)

base_model$md<-as.numeric(base_model$md)
base_model$only_provider<-as.numeric(base_model$only_provider)

base_model$personal_income<-as.numeric(base_model$personal_income)


mediate(emotional_e+depers+personal_a~ 
          md+
          work_hours+
          only_provider+ personal_income, data = base_model, n.iter = 10000) %>%
  print(large)


Mediation/Moderation Analysis 
Call: mediate(y = emotional_e + depers + personal_a ~ md + work_hours + 
    only_provider + personal_income, data = base_model, n.iter = 10000)

The DV (Y) was  emotional_e depers personal_a . The IV (X) was  md work_hours only_provider personal_income . The mediating variable(s) =  .Call: mediate(y = emotional_e + depers + personal_a ~ md + work_hours + 
    only_provider + personal_income, data = base_model, n.iter = 10000)

No mediator specified leads to traditional regression 
                emotional_e   se     t  df     Prob
Intercept              0.32 0.04  8.03 610 5.04e-15
md                     0.15 0.02  7.04 610 5.22e-12
work_hours             0.00 0.00  3.79 610 1.66e-04
only_provider         -0.02 0.02 -0.95 610 3.43e-01
personal_income       -0.04 0.01 -3.14 610 1.75e-03

R = 0.33 R2 = 0.11   F = 19.23 on 4 and 610 DF   p-value:  6.57e-15 
                depers   se     t  df     Prob
Intercept         0.20 0.03  7.18 610 2.05e-12
md                0.05 0.02  3.17 610 1.62e-03
work_hours        0.00 0.00  2.05 610 4.08e-02
only_provider     0.01 0.01  0.61 610 5.44e-01
personal_income  -0.01 0.01 -0.84 610 4.04e-01

R = 0.17 R2 = 0.03   F = 4.8 on 4 and 610 DF   p-value:  0.000803 
                personal_a   se     t  df     Prob
Intercept             0.67 0.03 22.25 610 9.48e-81
md                   -0.02 0.02 -1.30 610 1.94e-01
work_hours            0.00 0.00 -0.85 610 3.93e-01
only_provider         0.03 0.02  2.24 610 2.56e-02
personal_income       0.01 0.01  1.40 610 1.62e-01

R = 0.13 R2 = 0.02   F = 2.77 on 4 and 610 DF   p-value:  0.0264