Introducción general

Objetivo

Comprobar los supuestos de la logoterapia de que el Sindrome de Burnout está asociado a una ausencia de sentido de la vida y un proyecto de vida laboral inauténtico.

Fuente

Cuestionario autoadministrado, anónimo y de respuesta voluntaria (estudio transversal).

Población bajo estudio

Médicos/as residentes de algunos programas del Hospital Italiano de Buenos Aires.

Tamaño de la muestra

104, sin datos faltantes.

Variable respuesta

Presencia de sindrome de burnout (BO), dicotómica, según Maslach Burnout Inventory; si 3 escalas positivas: escala de agotamiento emocional (BO si >21), escala de despersonalización (BO si >6) y escala de realización personal (BO si <36).

Variables explicativas

Proyect in Life Test (PIL test; variable continua): mide “sentido de la vida”.

Falta de sentido de la vida (“vacío existencial”; variable dicotómica): si PIL test <105.

Cuestionario de proyecto de vida laboral (variable continua): a mayor valor, mayor autenticidad del proyecto de vida laboral.

Proyecto de vida laboral inauténtico (variable dicotómica): si cuestionario de proyecto de vida laboral <36.

Sexo (variable dicotómica)

Edad (variable continua)

Año de residencia (variable ordinal)

Especialidad (variable categórica politómica)

Especialidad clínica (variable dicotómica)

Carga de datos y de librerías

Primero cargamos las librerías necesarias:

library(tidyverse)
library(dplyr)
library(readxl)
library(car)
library(GGally)
library(arsenal)
library(kableExtra)
library(caret)
library(coin)
library(glmnet)
library(ROCit)
library(ggeffects)
library(ggpubr)
library(DHARMa)
library(gtsummary)
library(sjPlot)

Cargamos la base de datos original, modificando los nombres y eligiendo las variables de interés para la investigación:

datos <- read_excel("~/Downloads/BURNOUT/Logoterapia_base.xlsx")
datos$edad <- datos$Edad
datos$sexo <- datos$Sexo
datos$año <- datos$`Año de residencia`
datos$añonum <- as.numeric(datos$año)
datos$agot <- datos$`AGOTAMIENTO EMOCIONAL (SUMA TOTAL)`
datos$agot_dicot <- datos$`AGOTAMIENTO EMOCIONAL (NEIRA) MAYOR A 21 (SI O NO)`
datos$desp <- datos$`DESPERSONALIZACIÓN (SUMA TOTAL)`
datos$desp_dicot <- datos$`DESPERSONALIZACIÓN (NEIRA) MAYOR A 6 (SI O NO)`
datos$real <- datos$`REALIZACIÓN PERSONAL (SUMA TOTAL)`
datos$real_dicot <- datos$`ALTERACIÓN EN LA REALIZACIÓN PERSONAL (NEIRA) MENOR A 36  (SI O NO)`
datos$bo <- datos$`PRESENCIA DE BURNOUT 3/3 (NEIRA) ESCALAS POSITIVAS (SI O NO)`
datos$bo2 <- datos$`PRESENCIA DE BURNOUT 2/3 (NEIRA) ESCALAS POSITIVAS (SI O NO)`
datos$pil <- datos$`PIL TEST (SUMA TOTAL)`
datos$sentido_dicot <- datos$`FALTA DE UN CLARO SENTIDO O ZONA DE INDEFINICIÓN RESPECTO AL SENTIDO DE LA VIDA (MENOR A 106) (SI O NO)`
datos$proy <- datos$`PROYECTO DE VIDA LABORAL (SUMA TOTAL)`
datos$proy_dicot <- datos$`PROYECTO DE VIDA LABORAL INAUTÉNTICO (MENOR A 36) (SI O NO)`
datos$sobreadapt <- datos$`AFIRMATIVAS DE LA 1 A LA 5 (SOBREADAPTACIÓN, TOTAL)`
datos$inadapt <- datos$`AFIRMATIVAS DE LA 6 A LA 10 (INADAPTACIÓN, TOTAL)`
datos$adapt_dicot <- datos$`PREDOMINANTEMENTE CON SOBREADAPTACIÓN O INADAPTACIÓN`
datos$especialidad <- datos$Especialidad
datos$clinica<-0
datos$clinica[datos$Especialidad=="Clínica"]<-1
datos <- datos[,88:108]
datos$año <- as.factor(datos$año)
datos$sexo <- as.factor(datos$sexo)
datos$agot_dicot <- as.factor(datos$agot_dicot)
datos$desp_dicot <- as.factor(datos$desp_dicot)
datos$real_dicot <- as.factor(datos$real_dicot)
datos$bo <- as.factor(datos$bo)
datos$bo2 <- as.factor(datos$bo2)
datos$sentido_dicot <- as.factor(datos$sentido_dicot)
datos$proy_dicot <- as.factor(datos$proy_dicot)
datos$especialidad <- as.factor(datos$especialidad)
datos$clinica <- as.factor(datos$clinica)

labels(datos) <- c(edad="Edad (años)", sexo = "Sexo", año="Año de residencia",
                   agot="Agotamiento (total)", desp="Despersonalización (total)",
                   real="Realización personal (total)", pil="PIL test (total)",
                   proy="Proyecto de vida laboral (total)", clinica= "Especialidad clínica", 
                   especialidad= "Especialidad",bo="Presencia de burnout",
                   sentido_dicot="Falta de sentido de la vida", proy_dicot="Proyecto de vida laboral inauténtico")
levels(datos$sexo) <- c("Femenino", "Masculino")
levels(datos$bo) <- c("No", "Si")
levels(datos$clinica)<- c("No", "Si")
levels(datos$sentido_dicot)<- c("No", "Si")
levels(datos$proy_dicot) <- c("No", "Si")

Análisis exploratorio de variables

Burnout

30 de 104 encuestados (28.8%) presentaron diagnóstico de BO:

kable(table(datos$bo), col.names = c("Burnout","Total"), align ="l")%>%kable_styling()
Burnout Total
No 74
Si 30 
kable(round(prop.table(table(datos$bo))*100,2), col.names = c("Burnout","Porcentaje"),align="l")%>%kable_styling()
Burnout Porcentaje
No 71.15
Si 28.85

Edad:

La distribución de los residuos es no normal y homocedástica, con un boxplot con distribuciones similares entre los grupos con y sin burnout, por lo que usamos el test de Wilcoxon para comparar grupos.

modedad <- lm(edad~bo, datos)
qqnorm(modedad$residuals, pch=20)
qqline(modedad$residuals)

shapiro.test(modedad$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  modedad$residuals
## W = 0.9422, p-value = 0.0001947
leveneTest(edad~bo, datos)
ggplot(datos, aes(bo, edad))+geom_boxplot(show.legend=F, fill=c("chartreuse", "firebrick"))+ ggtitle("Boxplot de presencia de burnout según edad")+xlab("Presencia de burnout")+ylab("Edad (años)")+ theme_pubr()

wilcox.test(edad~bo, datos)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  edad by bo
## W = 1144, p-value = 0.8078
## alternative hypothesis: true location shift is not equal to 0

Sexo:

Usamos test de chi-cuadrado por presentar valores esperados mayores a 5 en todas las casillas.

chisq.test(datos$sexo, datos$bo)$expected
##            datos$bo
## datos$sexo        No       Si
##   Femenino  47.67308 19.32692
##   Masculino 26.32692 10.67308
chisq.test(datos$sexo, datos$bo)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  datos$sexo and datos$bo
## X-squared = 0.28127, df = 1, p-value = 0.5959

Año de residencia:

Usamos test de Fisher por no cumplir con los supuestos de valores esperados del test de chi cuadrado:

chisq.test(datos$año, datos$bo)$expected
##          datos$bo
## datos$año         No         Si
##         1 27.7500000 11.2500000
##         2 13.5192308  5.4807692
##         3 15.6538462  6.3461538
##         4 10.6730769  4.3269231
##         5  5.6923077  2.3076923
##         6  0.7115385  0.2884615
fisher.test(datos$año, datos$bo)
## 
##  Fisher's Exact Test for Count Data
## 
## data:  datos$año and datos$bo
## p-value = 0.8374
## alternative hypothesis: two.sided

Especialidad:

También usamos test de Fisher: la diferencia entre grupos es significativa. Vemos que la distribución de los grupos es muy desigual con algunas categorías vacias.

kable(table(datos$especialidad,datos$bo), col.names = c("Burnout No","Burnout Sí"), align ="l")%>%kable_styling()
Burnout No Burnout Sí
Cardiología 10 5
Clínica 33 23
Ginecología 5 0
Medicina Familiar 3 2
Pediatría 15 0
Terapia intensiva 8 0 
chisq.test(datos$especialidad, datos$bo)$expected
##                    datos$bo
## datos$especialidad         No        Si
##   Cardiología       10.673077  4.326923
##   Clínica           39.846154 16.153846
##   Ginecología        3.557692  1.442308
##   Medicina Familiar  3.557692  1.442308
##   Pediatría         10.673077  4.326923
##   Terapia intensiva  5.692308  2.307692
fisher.test(datos$especialidad, datos$bo)
## 
##  Fisher's Exact Test for Count Data
## 
## data:  datos$especialidad and datos$bo
## p-value = 0.002341
## alternative hypothesis: two.sided

Al dicotomizar en clínica/no clínica la diferencia también es significativa:

kable(table(datos$clinica,datos$bo), col.names = c("Burnout No","Burnout Sí"), row.names=T,align ="l")%>%kable_styling()
Burnout No Burnout Sí
No 41 7
Si 33 23 
chisq.test(datos$clinica, datos$bo)$expected
##              datos$bo
## datos$clinica       No       Si
##            No 34.15385 13.84615
##            Si 39.84615 16.15385
chisq.test(datos$clinica, datos$bo)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  datos$clinica and datos$bo
## X-squared = 7.5917, df = 1, p-value = 0.005864

PIL (Project in life test):

Distribución no normal homocedástica, con boxplots por grupo muy distintos, usamos test de la mediana, con diferencia significativa entre grupos.

modpil <- lm(pil~bo, datos)
qqnorm(modpil$residuals)
qqline(modpil$residuals)

shapiro.test(modpil$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  modpil$residuals
## W = 0.94794, p-value = 0.0004591
leveneTest(pil~bo, datos)
ggplot(datos, aes(bo, pil))+geom_boxplot(show.legend=F, fill=c("chartreuse", "firebrick"))+ ggtitle("Boxplot de presencia de burnout según puntaje total de PIL test")+xlab("Presencia de burnout")+ylab("PIL test total")+ theme_pubr()

median_test(pil~bo, datos)
## 
##  Asymptotic Two-Sample Brown-Mood Median Test
## 
## data:  pil by bo (No, Si)
## Z = 3.6313, p-value = 0.000282
## alternative hypothesis: true mu is not equal to 0

Ausencia de sentido de la vida (dicotómica)

El test de chi cuadrado es significativo al comparar grupos sin y con burnout.

chisq.test(datos$sentido_dicot, datos$bo)$expected
##                    datos$bo
## datos$sentido_dicot       No        Si
##                  No 60.48077 24.519231
##                  Si 13.51923  5.480769
chisq.test(datos$sentido_dicot, datos$bo)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  datos$sentido_dicot and datos$bo
## X-squared = 25.522, df = 1, p-value = 4.374e-07

Proyecto de vida laboral (puntaje total del cuestionario):

Distribución no normal homocedástica, con boxplots por grupo muy distintos, usamos test de la mediana, con diferencia significativa entre grupos.

modproy <- lm(proy~bo, datos)
qqnorm(modproy$residuals)
qqline(modproy$residuals)

shapiro.test(modproy$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  modproy$residuals
## W = 0.96459, p-value = 0.006985
#Proyecto no normal, uso wilcoxon
ggplot(datos, aes(bo, proy))+geom_boxplot(show.legend=F, fill=c("chartreuse", "firebrick"))+ ggtitle("Boxplot de presencia de burnout según puntaje total de test de proyecto de vida laboral")+xlab("Presencia de burnout")+ylab("Proyecto de vida laboral (total del test)")+ theme_pubr()

median_test(proy~bo, datos)
## 
##  Asymptotic Two-Sample Brown-Mood Median Test
## 
## data:  proy by bo (No, Si)
## Z = 4.308, p-value = 1.648e-05
## alternative hypothesis: true mu is not equal to 0

Proyecto de vida laboral inauténtico (dicotómica):

El test de chi cuadrado es significativo al comparar grupos sin y con burnout.

chisq.test(datos$bo, datos$proy_dicot)$expected
##         datos$proy_dicot
## datos$bo       No        Si
##       No 51.23077 22.769231
##       Si 20.76923  9.230769
chisq.test(datos$bo, datos$proy_dicot)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  datos$bo and datos$proy_dicot
## X-squared = 33.105, df = 1, p-value = 8.729e-09

Tablas

Plasmamos los resultados previos en una Tabla 1:

tabla1 <- tableby(bo~wt(edad, "median", "iqr")+chisq(sexo)+
                    fe(año)+fe(especialidad)+chisq(clinica)+medtest(pil, "median", "iqr")+ chisq(sentido_dicot)+medtest(proy, "median", "iqr")+chisq(proy_dicot), data=datos,
                  control=tableby.control(
                    stats.labels = list(median="Mediana", iqr="IIC"), test.pname ="p-valor",digits=2))
summary(tabla1, text=T, pfootnote = T) %>%
  kbl(caption = "Tabla 1: Variables explicativas según presencia o no de burnout")%>%
  row_spec(c(14,21,24,27,30,33), bold=T)%>%
  footnote(
    number=c("Test de Wilcoxon", "Test de Chi Cuadrado", "Test de Fisher", "Test de la mediana"))
Tabla 1: Variables explicativas según presencia o no de burnout
No (N=74) Si (N=30) Total (N=104) p-valor
Edad (años) 0.805 (1)
  • Mediana
27.00 27.50 27.00
  • IIC
3.00 4.00 4.00
Sexo 0.449 (2)
  • Femenino
46 (62.2%) 21 (70.0%) 67 (64.4%)
  • Masculino
28 (37.8%) 9 (30.0%) 37 (35.6%)
Año de residencia 0.837 (3)
  • 1
28 (37.8%) 11 (36.7%) 39 (37.5%)
  • 2
13 (17.6%) 6 (20.0%) 19 (18.3%)
  • 3
16 (21.6%) 6 (20.0%) 22 (21.2%)
  • 4
9 (12.2%) 6 (20.0%) 15 (14.4%)
  • 5
7 (9.5%) 1 (3.3%) 8 (7.7%)
  • 6
1 (1.4%) 0 (0.0%) 1 (1.0%)
Especialidad 0.002 (3)
  • Cardiología
10 (13.5%) 5 (16.7%) 15 (14.4%)
  • Clínica
33 (44.6%) 23 (76.7%) 56 (53.8%)
  • Ginecología
5 (6.8%) 0 (0.0%) 5 (4.8%)
  • Medicina Familiar
3 (4.1%) 2 (6.7%) 5 (4.8%)
  • Pediatría
15 (20.3%) 0 (0.0%) 15 (14.4%)
  • Terapia intensiva
8 (10.8%) 0 (0.0%) 8 (7.7%)
Especialidad clínica 0.003 (2)
  • No
41 (55.4%) 7 (23.3%) 48 (46.2%)
  • Si
33 (44.6%) 23 (76.7%) 56 (53.8%)
PIL test (total) < 0.001 (4)
  • Mediana
119.00 105.50 116.00
  • IIC
11.75 17.50 14.25
Falta de sentido de la vida < 0.001 (2)
  • No
70 (94.6%) 15 (50.0%) 85 (81.7%)
  • Si
4 (5.4%) 15 (50.0%) 19 (18.3%)
Proyecto de vida laboral (total) < 0.001 (4)
  • Mediana
37.00 34.00 36.50
  • IIC
4.00 3.75 4.00
Proyecto de vida laboral inauténtico < 0.001 (2)
  • No
64 (86.5%) 8 (26.7%) 72 (69.2%)
  • Si
10 (13.5%) 22 (73.3%) 32 (30.8%)
1 Test de Wilcoxon
2 Test de Chi Cuadrado
3 Test de Fisher
4 Test de la mediana

Confecciamos una Tabla 2 con los odds ratios de burnout crudos (resultantes de modelos logísticos simples) para cada variable, donde vemos que las variables significativas son especialidad clínica, PIL test (y su dicotómica) y Proyecto de vida laboral (y su dicotómica):

tbl_uv <-
  tbl_uvregression(
    datos %>% dplyr::select(bo, edad, sexo, año, clinica, pil, sentido_dicot, proy, proy_dicot),
    method = glm,
    y = bo,
    method.args = list(family = binomial),
    exponentiate = TRUE,
hide_n = TRUE)

tbl_uv %>% bold_labels()%>%
bold_p()%>%modify_header(label="**Variable**", ci="**IC 95%**", p.value="**p-valor**")%>% modify_caption("Tabla 2: Odds ratios crudos (según modelos logísticos simples)")
Tabla 2: Odds ratios crudos (según modelos logísticos simples)
Variable OR1 IC 95%1 p-valor
Edad (años) 0.96 0.79, 1.16 0.7
Sexo
Femenino — —
Masculino 0.70 0.27, 1.72 0.5
Año de residencia
1 — —
2 1.17 0.34, 3.83 0.8
3 0.95 0.28, 3.02 >0.9
4 1.70 0.47, 5.91 0.4
5 0.36 0.02, 2.40 0.4
6 0.00 >0.9
Especialidad clínica
No — —
Si 4.08 1.62, 11.4 0.004
PIL test (total) 0.90 0.85, 0.94 <0.001
Falta de sentido de la vida
No — —
Si 17.5 5.50, 68.5 <0.001
Proyecto de vida laboral (total) 0.64 0.51, 0.77 <0.001
Proyecto de vida laboral inauténtico
No — —
Si 17.6 6.45, 53.3 <0.001
1 OR = Odds Ratio, CI = Confidence Interval

Modelos de regresión

Construiremos diferentes modelos logísticos múltiples que expliquen a la variable de respuesta burnout. Considerando la “regla del pulgar”, la cantidad de variables explicativas que podamos sumar al modelo va a ser limitada (1 variable explicativa por cada 10 eventos). En la muestra de 104 residentes evaluados 30 presentaron el evento, por lo tanto, nuestro modelo podrá tener hasta 3 variables explicativas.

Modelo de 3 variables

Empezamos entonces armando un modelo logístico múltiple incluyendo las 3 variables explicativas dicotómicas (especialidad clínica, ausencia de sentido de la vida, proyecto laboral inauténtico) que resultaron significativamente asociadas a la presencia de burnout en los modelos simples:

modclinica <- glm(bo~clinica+sentido_dicot+proy_dicot,family=binomial,datos)
summary(modclinica)
## 
## Call:
## glm(formula = bo ~ clinica + sentido_dicot + proy_dicot, family = binomial, 
##     data = datos)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.1021  -0.5607  -0.3609   0.4822   2.3510  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      -2.6984     0.5328  -5.065 4.08e-07 ***
## clinicaSi         0.9277     0.5903   1.571 0.116071    
## sentido_dicotSi   1.6862     0.7529   2.240 0.025111 *  
## proy_dicotSi      2.1777     0.5958   3.655 0.000257 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 124.96  on 103  degrees of freedom
## Residual deviance:  80.90  on 100  degrees of freedom
## AIC: 88.9
## 
## Number of Fisher Scoring iterations: 5
exp(modclinica$coefficients)
##     (Intercept)       clinicaSi sentido_dicotSi    proy_dicotSi 
##      0.06731227      2.52868953      5.39905477      8.82576201
exp(confint(modclinica))
##                     2.5 %     97.5 %
## (Intercept)     0.0209224  0.1729989
## clinicaSi       0.8090615  8.4667748
## sentido_dicotSi 1.2660690 25.5954989
## proy_dicotSi    2.7921483 29.5879865
plot_model(modclinica, show.values = T,auto.label = T, prefix.labels = "label")+ggtitle("Odds ratios para presencia de burnout", subtitle="En modelo de 3 variables")+theme_pubr()

Vemos que la variable “especialidad clínica”, al ser controlada por las otras dos variables, deja de ser significativa. Recordemos la distribución notoriamente desigual entre especialidades: puede ser que nuestra muestra esté sesgada en cuanto a la variable “clínica”, por lo tanto no la tendremos en cuenta para los modelos finales. Las variables que restan son los totales numéricos de ambos tests, o las variables categóricas correspondientes (ausencia de sentido de la vida, y proyecto de vida laboral inauténtico). Probemos primero el modelo con ambas variables categóricas.

Modelo de 2 variables categóricas

mod <- glm(bo~sentido_dicot+proy_dicot, family=binomial, datos)
summary(mod)
## 
## Call:
## glm(formula = bo ~ sentido_dicot + proy_dicot, family = binomial, 
##     data = datos)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.9942  -0.4558  -0.4558   0.5427   2.1521  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      -2.2119     0.3893  -5.682 1.33e-08 ***
## sentido_dicotSi   1.7935     0.7205   2.489 0.012798 *  
## proy_dicotSi      2.2596     0.5831   3.875 0.000106 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 124.960  on 103  degrees of freedom
## Residual deviance:  83.442  on 101  degrees of freedom
## AIC: 89.442
## 
## Number of Fisher Scoring iterations: 4

Vemos que ambas variables resultan significativas; veamos si su interacción también lo es:

modint <- glm(bo~sentido_dicot*proy_dicot, family=binomial, datos)
Anova(modint)

Como la interacción entre ambas variables categóricas no es significativa, tiene sentido informar los OR de efectos simples:

exp(mod$coefficients)
##     (Intercept) sentido_dicotSi    proy_dicotSi 
##       0.1094913       6.0102355       9.5791221
exp(confint(mod))
##                      2.5 %   97.5 %
## (Intercept)     0.04694411  0.22019
## sentido_dicotSi 1.51532836 26.90209
## proy_dicotSi    3.11228174 31.31051
tbl_dicot <- tbl_regression(mod, exponentiate = T)

tbl_dicot %>% bold_labels()%>%
bold_p()%>%modify_header(label="**Variable**", ci="**IC 95%**", p.value="**p-valor**")%>% modify_caption("Tabla 3: Odds ratios ajustados (según modelo con 2 variables dicotómicas)")
Tabla 3: Odds ratios ajustados (según modelo con 2 variables dicotómicas)
Variable OR1 IC 95%1 p-valor
Falta de sentido de la vida
No — —
Si 6.01 1.52, 26.9 0.013
Proyecto de vida laboral inauténtico
No — —
Si 9.58 3.11, 31.3 <0.001
1 OR = Odds Ratio, CI = Confidence Interval 
plot_model(mod, show.values = T,auto.label = T, prefix.labels = "label")+ggtitle("Odds ratios para presencia de burnout", subtitle="En modelo de 2 variables dicotómicas")+theme_pubr()

Vemos una gran magnitud de efecto de ambas variables explicativas: los residentes con ausencia del sentido de la vida tienen un odds de tener burnout de 6 veces el odds de los que tienen valores normales del PIL test (IC 1.5-26.9); y quienes tienen un proyecto de vida laboral inauténtico tienen un odds de burnout de 9.6 veces en relación a quienes tienen un proyecto de vida laboral auténtico (IC 3.1 - 31.3). Se ve también que os intervalos de confianza son muy grandes.

Chequeamos los supuestos del modelo mediante la librería DHARMa (residual diagnostics for hierarchical (multi-level/mixed) regression models):

simulationOutput <- simulateResiduals(fittedModel = mod, plot = T)

Los supuestos del modelo se cumplen adecuadamente.

Confeccionamos la curva ROC de éste modelo mediante validación cruzada múltiple (5-fold, 500 repeticiones), calculando su AUC; y calculamos la sensibilidad y especificidad en el punto óptimo de trade-off entre ambas:

set.seed(123)
ctrl <- trainControl(method="repeatedcv", number=5, repeats=500, classProbs = TRUE)
modcaret <- train(bo~sentido_dicot+proy_dicot, data=datos,method="glm", trControl=ctrl)
print(modcaret)
## Generalized Linear Model 
## 
## 104 samples
##   2 predictor
##   2 classes: 'No', 'Si' 
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold, repeated 500 times) 
## Summary of sample sizes: 83, 83, 83, 84, 83, 84, ... 
## Resampling results:
## 
##   Accuracy   Kappa   
##   0.7945724  0.473817
summary(modcaret)
## 
## Call:
## NULL
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.9942  -0.4558  -0.4558   0.5427   2.1521  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      -2.2119     0.3893  -5.682 1.33e-08 ***
## sentido_dicotSi   1.7935     0.7205   2.489 0.012798 *  
## proy_dicotSi      2.2596     0.5831   3.875 0.000106 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 124.960  on 103  degrees of freedom
## Residual deviance:  83.442  on 101  degrees of freedom
## AIC: 89.442
## 
## Number of Fisher Scoring iterations: 4
ajustadoscaret<-predict(modcaret,newdata = datos,type="prob")[2]
ROC <- rocit(score=ajustadoscaret$Si,class=ifelse(datos$bo=="Si",1,0))
plot(ROC, col=3, xlab=" Tasa de falsos positivos (1-Especificidad)", 
     ylab="Sensibilidad (Tasa de verdaderos positivos)")

ROC$AUC
## [1] 0.8315315
distancia <- sqrt(ROC$FPR^2+(ROC$TPR-1)^2)
minimo <- ROC$Cutoff[which.min(distancia)]
predichoscaret<-as.factor(ifelse(ajustadoscaret$Si<minimo,"No","Si"))
tablacaret<-table(Observados=datos$bo, Predichos=predichoscaret)
tablacaret
##           Predichos
## Observados No Si
##         No 62 12
##         Si  7 23
sensibilidad <- tablacaret[2,2]/(tablacaret[2,2]+tablacaret[2,1])
especificidad <- tablacaret[1,1]/(tablacaret[1,1]+tablacaret[1,2])
sensibilidad
## [1] 0.7666667
especificidad
## [1] 0.8378378

Entonces, el modelo con las 2 variables explicativas dicotómicas tiene un AUC de 0.83, una sensibilidad de 76.7% y una especificidad del 83.8%.

Modelo de 2 variables continuas

Probemos ahora el modelo con ambas variables continuas:

mod2 <- glm(bo~pil+proy, family=binomial, datos)
summary(mod2)
## 
## Call:
## glm(formula = bo ~ pil + proy, family = binomial, data = datos)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.6035  -0.6902  -0.3982   0.5048   2.0216  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 16.65375    3.87379   4.299 1.72e-05 ***
## pil         -0.05376    0.03218  -1.671   0.0947 .  
## proy        -0.32139    0.11625  -2.765   0.0057 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 124.960  on 103  degrees of freedom
## Residual deviance:  88.305  on 101  degrees of freedom
## AIC: 94.305
## 
## Number of Fisher Scoring iterations: 5

Vemos que una de ambas variables (puntaje total del PIL test) no resulta significativa. Analicemos la correlación entre ambas variables:

cor(datos$pil, datos$proy)
## [1] 0.6711153
ggplot(datos, aes(pil, proy))+geom_point(col="chocolate")+xlab("Puntaje total de PIL test")+ylab("Puntaje total de cuestionario de proyecto de vida laboral")

Vemos que la correlación entre ambas variables es grande. Veamos que de todas formas el VIF del modelo no es preocupante:

library(car)
vif(mod2)
##      pil     proy 
## 1.347787 1.347787

Para evaluar la posibilidad de interacción entre ambas variables explicativas numéricas, confeccionamos gráficos de perfiles (previamente construyendo 3 categorías en cada variable):

quantile(datos$pil, seq(0,1,0.33333))
##       0%  33.333%  66.666%  99.999% 
##  74.0000 111.3330 120.0000 134.9969
datos$pil_categ <- factor(ifelse(datos$pil<111.33, "bajo", ifelse(
  datos$pil<120, "medio", "alto"
)),levels=c("bajo", "medio", "alto"))

quantile(datos$proy, seq(0,1,0.33333))
##       0%  33.333%  66.666%  99.999% 
## 25.00000 36.00000 37.00000 44.99897
datos$proy_categ <-factor(ifelse(datos$proy<36, "bajo", ifelse(
  datos$proy<38, "medio", "alto"
)), levels=c("bajo", "medio", "alto"))

medias <- aggregate(pil~proy_categ+bo, datos, mean)
ggplot(medias, aes(x=proy_categ, y=pil, group=bo, col=bo))+
  geom_line(aes(linetype=bo))+
  geom_point(aes(shape=bo))+
  theme_pubr()

ggplot(medias, aes(x=bo, y=pil, group=proy_categ, col=proy_categ))+
  geom_line(aes(linetype=proy_categ))+
  geom_point(aes(shape=proy_categ))+
  theme_pubr()

medias <- aggregate(proy~pil_categ+bo, datos, mean)
ggplot(medias, aes(x=pil_categ, y=proy, group=bo, col=bo))+
  geom_line(aes(linetype=bo))+
  geom_point(aes(shape=bo))+
  theme_pubr()

ggplot(medias, aes(x=bo, y=proy, group=pil_categ, col=pil_categ))+
  geom_line(aes(linetype=pil_categ))+
  geom_point(aes(shape=pil_categ))+
  theme_pubr()

Se observa que en algunos gráficos cambian las tendencias segun el valor de la otra variable, sugiriendo una interacción. Probemos entonces un modelo con interacción entre ambas variables:

mod3 <- glm(bo~pil*proy, family=binomial, datos)
summary(mod3)
## 
## Call:
## glm(formula = bo ~ pil * proy, family = binomial, data = datos)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.8901  -0.6803  -0.1951   0.6699   2.4031  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)   
## (Intercept) -68.11800   28.16769  -2.418  0.01559 * 
## pil           0.75481    0.28025   2.693  0.00707 **
## proy          2.21210    0.86258   2.565  0.01033 * 
## pil:proy     -0.02401    0.00842  -2.852  0.00435 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 124.960  on 103  degrees of freedom
## Residual deviance:  79.212  on 100  degrees of freedom
## AIC: 87.212
## 
## Number of Fisher Scoring iterations: 6
Anova(mod3)

Vemos que la interacción es significativa y desciende el AIC (de hecho es el menor AIC de los modelos vistos hasta el momento). Podemos visualizar el efecto de la interacción graficando la probabilidad de presentar burnout según cada variable, dados 3 valores de la otra variable (mediana +/- 1 desvío standard):

modgg <- glm(bo~pil*proy, family=binomial, datos)
interaccion1 <- ggpredict(modgg, terms=c("proy","pil"))
interaccion2<- ggpredict(modgg, terms=c("pil", "proy"))
a <- plot(interaccion1,add.data = T,jitter = 0.1, show.title = F, alpha = 0.1)
a+theme_pubr()+ylab("Probabilidad de presentar burnout")

b<- plot(interaccion2,add.data = T,jitter = 0.1, show.title = F, alpha=0.1)
b+theme_pubr()+ylab("Probabilidad de presentar burnout")

Se observa que a valores medios y altos de PIL test, presentar valores bajos de proyecto de vida laboral aumenta notoriamente la probabilidad de presentar burnout, y viceversa.

Chequeemos los supuestos de este modelo:

simulationOutput <- simulateResiduals(fittedModel = mod3, seed=123, plot = T)

Validemos entonces este modelo con curva ROC, AUC, sensibilidad y especificidad:

set.seed(123)
ctrl <- trainControl(method="repeatedcv", number=5, repeats=500, classProbs = TRUE)
modcaret <- train(bo~pil*proy, data=datos,method="glm", trControl=ctrl)
print(modcaret)
## Generalized Linear Model 
## 
## 104 samples
##   2 predictor
##   2 classes: 'No', 'Si' 
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold, repeated 500 times) 
## Summary of sample sizes: 83, 83, 83, 84, 83, 84, ... 
## Resampling results:
## 
##   Accuracy   Kappa    
##   0.8280981  0.5658208
summary(modcaret)
## 
## Call:
## NULL
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.8901  -0.6803  -0.1951   0.6699   2.4031  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)   
## (Intercept) -68.11800   28.16769  -2.418  0.01559 * 
## pil           0.75481    0.28025   2.693  0.00707 **
## proy          2.21210    0.86258   2.565  0.01033 * 
## `pil:proy`   -0.02401    0.00842  -2.852  0.00435 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 124.960  on 103  degrees of freedom
## Residual deviance:  79.212  on 100  degrees of freedom
## AIC: 87.212
## 
## Number of Fisher Scoring iterations: 6
ajustadoscaret<-predict(modcaret,newdata = datos,type="prob")[2]
ROC <- rocit(score=ajustadoscaret$Si,class=ifelse(datos$bo=="Si",1,0))
plot(ROC, col=3, xlab=" Tasa de falsos positivos (1-Especificidad)", 
     ylab="Sensibilidad (Tasa de verdaderos positivos)")

ROC$AUC
## [1] 0.8752252
distancia <- sqrt(ROC$FPR^2+(ROC$TPR-1)^2)
minimo <- ROC$Cutoff[which.min(distancia)]
predichoscaret<-as.factor(ifelse(ajustadoscaret$Si<minimo,"No","Si"))
tablacaret<-table(Observados=datos$bo, Predichos=predichoscaret)
tablacaret
##           Predichos
## Observados No Si
##         No 59 15
##         Si  6 24
sensibilidad <- tablacaret[2,2]/(tablacaret[2,2]+tablacaret[2,1])
especificidad <- tablacaret[1,1]/(tablacaret[1,1]+tablacaret[1,2])
sensibilidad
## [1] 0.8
especificidad
## [1] 0.7972973

Vemos entonces que cuando incluimos ambas variables explicativas numéricas continuas y con interacción entre ambas, aumenta solo ligeramente el valor predictivo del modelo: mejora la sensibilidad a 80% y el AUC a 0.87, y empeora la especificidad bajando a 79.7%.