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 %>% select(origin,age, gender, family_income, personal_income,only_provider, laboral, work_type, work_md) %>%
  tbl_summary(by=origin) %>% 
  add_p() %>%
  add_q() %>%
  modify_caption("**Table 1** Demographics features for center (N = {N})")

Table 1 Demographics features for center (N = 927)

Characteristic	Eva Peron, N = 1051	Fleni, N = 6151	Mendoza, N = 281	Other, N = 1791	p-value2	q-value3
age	45 (35, 55)	40 (33, 51)	41 (35, 47)	46 (38, 57)	<0.001	<0.001
Unknown	7	348	1	13
gender					0.039	0.039
Male	19 (18%)	190 (31%)	9 (32%)	57 (32%)
Female	84 (80%)	423 (69%)	19 (68%)	122 (68%)
Other	2 (1.9%)	2 (0.3%)	0 (0%)	0 (0%)
family_income
Less than $100.000	0 (0%)	5 (0.8%)	1 (3.6%)	1 (0.6%)
$100.001 to $ 250.000	23 (22%)	104 (17%)	1 (3.6%)	8 (4.5%)
$250.001 to 400.000	35 (33%)	177 (29%)	5 (18%)	36 (20%)
More than $400.000	47 (45%)	329 (53%)	21 (75%)	134 (75%)
personal_income
Less than $100.000	1 (1.0%)	9 (1.5%)	0 (0%)	3 (1.7%)
$100.001 to $ 250.000	35 (33%)	196 (32%)	5 (18%)	23 (13%)
$250.001 to 400.000	39 (37%)	225 (37%)	8 (29%)	41 (23%)
More than $400.000	30 (29%)	185 (30%)	15 (54%)	112 (63%)
only_provider	72 (69%)	396 (64%)	19 (68%)	137 (77%)	0.025	0.033
laboral
Private practice	4 (3.8%)	519 (85%)	10 (36%)	93 (52%)
Public	62 (60%)	15 (2.4%)	3 (11%)	42 (23%)
Both mostly public	29 (28%)	18 (2.9%)	6 (21%)	19 (11%)
Both mostly private	9 (8.7%)	62 (10%)	9 (32%)	25 (14%)
Unknown	1	1	0	0
work_type
Public-facing customer service (receptionists and secretaries)	17 (16%)	78 (13%)	0 (0%)	6 (3.4%)
Direct patient care (physician, nurse, physiotherapist, therapist, etc.)	67 (64%)	369 (60%)	26 (93%)	162 (91%)
Indirect patient care (biochemist, imaging studies report, etc.)	15 (14%)	71 (12%)	2 (7.1%)	4 (2.2%)
Administrative tasks related to healthcare without public-facing duties (commercial administrative sector)	6 (5.7%)	97 (16%)	0 (0%)	7 (3.9%)
work_md	46 (44%)	189 (31%)	16 (57%)	142 (79%)	<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

Mucha gente no quizo dar su año de nacimiento por temor a ser identificados en Fleni, hay un numero muy grande de perdidos, no vamos a poder usar ese dato

base %>% select(origin,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 = 927)

Characteristic	No, N = 5341	Yes, N = 3931	p-value2	q-value3
origin			<0.001	<0.001
Eva Peron	59 (11%)	46 (12%)
Fleni	426 (80%)	189 (48%)
Mendoza	12 (2.2%)	16 (4.1%)
Other	37 (6.9%)	142 (36%)
age	42 (34, 51)	45 (36, 55)	0.002	0.002
Unknown	250	119
gender			<0.001	<0.001
Male	122 (23%)	153 (39%)
Female	409 (77%)	239 (61%)
Other	3 (0.6%)	1 (0.3%)
family_income			<0.001	<0.001
Less than $100.000	7 (1.3%)	0 (0%)
$100.001 to $ 250.000	124 (23%)	12 (3.1%)
$250.001 to 400.000	169 (32%)	84 (21%)
More than $400.000	234 (44%)	297 (76%)
personal_income			<0.001	<0.001
Less than $100.000	11 (2.1%)	2 (0.5%)
$100.001 to $ 250.000	228 (43%)	31 (7.9%)
$250.001 to 400.000	195 (37%)	118 (30%)
More than $400.000	100 (19%)	242 (62%)
only_provider	312 (58%)	312 (79%)	<0.001	<0.001
laboral			<0.001	<0.001
Private practice	392 (73%)	234 (60%)
Public	78 (15%)	44 (11%)
Both mostly public	21 (3.9%)	51 (13%)
Both mostly private	43 (8.1%)	62 (16%)
Unknown	0	2
work_type			<0.001	<0.001
Public-facing customer service (receptionists and secretaries)	100 (19%)	1 (0.3%)
Direct patient care (physician, nurse, physiotherapist, therapist, etc.)	272 (51%)	352 (90%)
Indirect patient care (biochemist, imaging studies report, etc.)	61 (11%)	31 (7.9%)
Administrative tasks related to healthcare without public-facing duties (commercial administrative sector)	101 (19%)	9 (2.3%)
1 n (%); Median (IQR)
2 Pearson’s Chi-squared test; Wilcoxon rank sum test; Fisher’s exact test
3 False discovery rate correction for multiple testing

base %>% 
  filter(work_md=="Yes") %>%
  select(origin,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 = 391)

Characteristic	Overall, N = 3911	Private, N = 2961	Public, N = 951	p-value2	q-value3
origin				<0.001	<0.001
Eva Peron	45 (12%)	10 (3.4%)	35 (37%)
Fleni	188 (48%)	179 (60%)	9 (9.5%)
Mendoza	16 (4.1%)	11 (3.7%)	5 (5.3%)
Other	142 (36%)	96 (32%)	46 (48%)
age	45 (36, 55)	44 (36, 55)	45 (36, 56)	0.9	0.9
Unknown	118	106	12
gender				0.5	0.7
Male	153 (39%)	120 (41%)	33 (35%)
Female	237 (61%)	175 (59%)	62 (65%)
Other	1 (0.3%)	1 (0.3%)	0 (0%)
years_of_graduation	16 (8, 30)	15 (8, 28)	19 (8, 31)	0.3	0.6
work_type				0.011	0.046
Public-facing customer service (receptionists and secretaries)	1 (0.3%)	1 (0.3%)	0 (0%)
Direct patient care (physician, nurse, physiotherapist, therapist, etc.)	351 (90%)	266 (90%)	85 (89%)
Indirect patient care (biochemist, imaging studies report, etc.)	30 (7.7%)	26 (8.8%)	4 (4.2%)
Administrative tasks related to healthcare without public-facing duties (commercial administrative sector)	9 (2.3%)	3 (1.0%)	6 (6.3%)
family_income				0.3	0.6
Less than $100.000	0 (0%)	0 (0%)	0 (0%)
$100.001 to $ 250.000	12 (3.1%)	8 (2.7%)	4 (4.2%)
$250.001 to 400.000	82 (21%)	58 (20%)	24 (25%)
More than $400.000	297 (76%)	230 (78%)	67 (71%)
personal_income				0.6	0.8
Less than $100.000	2 (0.5%)	1 (0.3%)	1 (1.1%)
$100.001 to $ 250.000	30 (7.7%)	21 (7.1%)	9 (9.5%)
$250.001 to 400.000	117 (30%)	89 (30%)	28 (29%)
More than $400.000	242 (62%)	185 (63%)	57 (60%)
only_provider	311 (80%)	235 (79%)	76 (80%)	0.9	0.9
laboral				<0.001	<0.001
Private practice	234 (60%)	234 (79%)	0 (0%)
Public	44 (11%)	0 (0%)	44 (46%)
Both mostly public	51 (13%)	0 (0%)	51 (54%)
Both mostly private	62 (16%)	62 (21%)	0 (0%)
work_hours	44 (36, 55)	45 (35, 55)	40 (36, 55)	0.7	0.8
work_number
0	2 (0.5%)	2 (0.7%)	0 (0%)
1	91 (23%)	79 (27%)	12 (13%)
2	115 (29%)	80 (27%)	35 (37%)
3	108 (28%)	80 (27%)	28 (29%)
4	45 (12%)	33 (11%)	12 (13%)
5	20 (5.1%)	18 (6.1%)	2 (2.1%)
6	7 (1.8%)	3 (1.0%)	4 (4.2%)
7	3 (0.8%)	1 (0.3%)	2 (2.1%)
work_teaching	164 (42%)	120 (41%)	44 (46%)	0.3	0.6
work_emergency	111 (28%)	80 (27%)	31 (33%)	0.3	0.6
work_uti	45 (12%)	32 (11%)	13 (14%)	0.4	0.7
1 n (%); Median (IQR)
2 Fisher’s exact test; Wilcoxon rank sum test; Pearson’s Chi-squared 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 = 927)

Characteristic	Overall, N = 9271	No, N = 5341	Yes, N = 3931	p-value2	q-value3
dass_stress	12 (8, 20)	12 (6, 20)	14 (8, 20)	0.11	0.3
Unknown	409	200	209
dass_stress_dx				0.6	0.8
Normal	171 (36%)	106 (35%)	65 (37%)
Mild	161 (34%)	108 (36%)	53 (30%)
Moderate	103 (22%)	63 (21%)	40 (23%)
Severe	28 (5.9%)	16 (5.4%)	12 (6.8%)
Extremely severe	12 (2.5%)	6 (2.0%)	6 (3.4%)
Unknown	452	235	217
dass_anxiety	4 (0, 10)	4 (0, 10)	4 (0, 8)	0.6	0.8
Unknown	249	131	118
dass_anxiety_dx				0.8	0.8
Normal	249 (52%)	141 (50%)	108 (56%)
Mild	47 (9.8%)	29 (10%)	18 (9.3%)
Moderate	82 (17%)	51 (18%)	31 (16%)
Severe	39 (8.2%)	24 (8.5%)	15 (7.7%)
Extremely severe	61 (13%)	39 (14%)	22 (11%)
Unknown	449	250	199
dass_depression	6 (2, 12)	6 (2, 12)	6 (2, 14)	0.049	0.3
Unknown	285	147	138
dass_depression_dx				0.6	0.8
Normal	258 (52%)	160 (55%)	98 (48%)
Mild	84 (17%)	47 (16%)	37 (18%)
Moderate	108 (22%)	59 (20%)	49 (24%)
Severe	28 (5.6%)	15 (5.1%)	13 (6.3%)
Extremely severe	20 (4.0%)	11 (3.8%)	9 (4.4%)
Unknown	429	242	187
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 = 927)

Characteristic	Overall, N = 9271	No, N = 5341	Yes, N = 3931	p-value2	q-value3
maslach_emotional_exhaustion	27 (17, 37)	24 (15, 34)	32 (22, 41)	<0.001	<0.001
Unknown	4	3	1
maslach_depersonalization	8.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)	29 (24, 34)	29 (24, 33)	0.6	0.6
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 = 927)

Characteristic	Overall, N = 9271	No, N = 5341	Yes, N = 3931	p-value2	q-value3
maslach_emotional_exhaustion_index	0.50 (0.31, 0.69)	0.44 (0.28, 0.63)	0.58 (0.41, 0.76)	<0.001	<0.001
Unknown	4	3	1
maslach_depersonalization_index	0.27 (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.69 (0.57, 0.81)	0.69 (0.57, 0.79)	0.6	0.6
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.28 0.03  9.93 923 3.88e-22
md                   0.11 0.02  6.73 923 2.97e-11
work_hours           0.00 0.00  5.02 923 6.08e-07
only_provider       -0.02 0.02 -1.32 923 1.86e-01

R = 0.29 R2 = 0.08   F = 27.76 on 3 and 923 DF   p-value:  3.43e-17 
              personal_a   se     t  df      Prob
Intercept           0.69 0.02 34.25 923 1.54e-166
md                 -0.01 0.01 -0.51 923  6.12e-01
work_hours          0.00 0.00 -0.55 923  5.80e-01
only_provider       0.02 0.01  1.80 923  7.24e-02

R = 0.06 R2 = 0   F = 1.16 on 3 and 923 DF   p-value:  0.323 
              depers   se     t  df     Prob
Intercept       0.19 0.02  9.73 923 2.39e-21
md              0.05 0.01  4.53 923 6.79e-06
work_hours      0.00 0.00  2.84 923 4.64e-03
only_provider   0.00 0.01 -0.02 923 9.84e-01

R = 0.19 R2 = 0.04   F = 11.63 on 3 and 923 DF   p-value:  1.74e-07 

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.14 0.02  6.03 920 2.38e-09
depress              0.30 0.04  6.94 920 7.30e-12
anx                  0.07 0.04  1.81 920 7.14e-02
stress               0.35 0.04  8.87 920 3.89e-18
md                   0.09 0.01  7.05 920 3.48e-12
work_hours           0.00 0.00  3.91 920 9.90e-05
only_provider       -0.02 0.01 -1.24 920 2.16e-01

R = 0.62 R2 = 0.39   F = 97.09 on 6 and 920 DF   p-value:  1.62e-94 
              personal_a   se     t  df      Prob
Intercept           0.72 0.02 35.88 920 4.61e-177
depress            -0.24 0.04 -6.73 920  3.08e-11
anx                 0.06 0.03  1.80 920  7.18e-02
stress             -0.03 0.03 -0.92 920  3.59e-01
md                  0.00 0.01  0.20 920  8.44e-01
work_hours          0.00 0.00 -0.06 920  9.51e-01
only_provider       0.02 0.01  1.52 920  1.30e-01

R = 0.29 R2 = 0.09   F = 14.28 on 6 and 920 DF   p-value:  1.32e-15 
              depers   se     t  df     Prob
Intercept       0.15 0.02  7.96 920 5.08e-15
depress         0.15 0.03  4.34 920 1.56e-05
anx             0.07 0.03  2.06 920 3.99e-02
stress          0.00 0.03 -0.15 920 8.78e-01
md              0.05 0.01  4.33 920 1.67e-05
work_hours      0.00 0.00  2.57 920 1.02e-02
only_provider   0.00 0.01  0.21 920 8.37e-01

R = 0.31 R2 = 0.09   F = 15.8 on 6 and 920 DF   p-value:  2.46e-17 

#  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.17 0.03  5.93 919 4.19e-09
depress                0.30 0.04  6.91 919 8.99e-12
anx                    0.06 0.04  1.57 919 1.16e-01
stress                 0.35 0.04  8.93 919 2.31e-18
md                     0.10 0.01  7.09 919 2.73e-12
work_hours             0.00 0.00  3.92 919 9.33e-05
only_provider         -0.01 0.01 -0.47 919 6.42e-01
personal_income       -0.02 0.01 -1.65 919 9.89e-02

R = 0.62 R2 = 0.39   F = 83.77 on 7 and 919 DF   p-value:  4.32e-94 
                depers   se     t  df     Prob
Intercept         0.16 0.02  6.87 919 1.21e-11
depress           0.15 0.03  4.33 919 1.66e-05
anx               0.07 0.03  1.98 919 4.82e-02
stress            0.00 0.03 -0.14 919 8.91e-01
md                0.05 0.01  4.11 919 4.36e-05
work_hours        0.00 0.00  2.58 919 1.02e-02
only_provider     0.00 0.01  0.37 919 7.08e-01
personal_income   0.00 0.01 -0.46 919 6.46e-01

R = 0.31 R2 = 0.09   F = 13.56 on 7 and 919 DF   p-value:  9.33e-17 

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.33 0.03  9.94 922 3.33e-22
md                     0.13 0.02  7.37 922 3.82e-13
work_hours             0.00 0.00  5.07 922 4.92e-07
only_provider          0.00 0.02 -0.04 922 9.72e-01
personal_income       -0.03 0.01 -2.95 922 3.29e-03

R = 0.3 R2 = 0.09   F = 23.16 on 4 and 922 DF   p-value:  2.9e-18 
                depers   se     t  df     Prob
Intercept         0.20 0.02  8.86 922 4.18e-18
md                0.06 0.01  4.66 922 3.58e-06
work_hours        0.00 0.00  2.85 922 4.49e-03
only_provider     0.01 0.01  0.52 922 6.06e-01
personal_income  -0.01 0.01 -1.33 922 1.83e-01

R = 0.2 R2 = 0.04   F = 9.17 on 4 and 922 DF   p-value:  2.86e-07 
                personal_a   se     t  df      Prob
Intercept             0.66 0.02 27.49 922 5.46e-122
md                   -0.02 0.01 -1.35 922  1.76e-01
work_hours            0.00 0.00 -0.57 922  5.70e-01
only_provider         0.01 0.01  0.81 922  4.16e-01
personal_income       0.02 0.01  2.09 922  3.68e-02

R = 0.09 R2 = 0.01   F = 1.97 on 4 and 922 DF   p-value:  0.0973