── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.4 ✔ readr 2.1.5
✔ forcats 1.0.0 ✔ stringr 1.5.1
✔ ggplot2 3.5.1 ✔ tibble 3.2.1
✔ lubridate 1.9.3 ✔ tidyr 1.3.1
✔ 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(haven)library(ggstatsplot)
You can cite this package as:
Patil, I. (2021). Visualizations with statistical details: The 'ggstatsplot' approach.
Journal of Open Source Software, 6(61), 3167, doi:10.21105/joss.03167
library(lme4)
Loading required package: Matrix
Attaching package: 'Matrix'
The following objects are masked from 'package:tidyr':
expand, pack, unpack
Hay varios actores que han sido evaluados por los sujetos de GRUPO=1. Cambiamos la estructura de los datos para adecuarla a formato tidy y poder hacer análisis de los datos posteriormente:
Los aciertos y errores de los sujetos se han recogido en variables que empiezan por TRI_ y que tienen un valor de 1 si el sujeto ha asignado el triaje correspondiente al actor y un 0 si no lo ha hecho. Vamos a cambiar la estructura de los datos para poder hacer análisis de los datos posteriormente:
No esperaba que hubiese valores fuera de rango: 4. Casi todos además han ocurrido en el GRUPO==1. Aparentemente el valor “4” se ha asignado cunado un sujeto no sabía qué decir y parecen haber ocurrido casi todos en el GRUPO 1. Parece que eso se merece su propia variable:
Desde luego los del GRUPO 1, dudan más que los del grupo 0. Me parece que podríais intentar explicarlo tal vez de este modo:
Había una vez un cienpiés que se movía con gracia y sin esfuerzo, asombrando a todos los animales con su coordinación perfecta. Un día, una rana, intrigada por la habilidad del centípedo, le preguntó: “¿Cómo logras coordinar todas tus patas sin tropezar nunca?” El centípedo, que nunca antes había pensado en cómo lo hacía, empezó a reflexionar sobre la pregunta. Intentó concentrarse en el movimiento de cada una de sus patas, pero al hacerlo, se confundió y tropezó, cayendo al suelo.
listaPorACTOR= dfTriaje %>%split(.$ACTOR) listaGraficos = listaPorACTOR%>%map2(names(listaPorACTOR), ~ {.x %>% ggstatsplot::ggbarstats(y=GRUPO, x=NOTRIAJE,bf.message =FALSE)+ggtitle(.y)})#Crear con la lista de Gráficos uno solo de dos columnaslistaGraficos
$ARTURO
$JAVIER
$LAURA
$MARIA
$MARIO
$MATEO
$NICOLAS
Vamos a estudiar la situación a nivel de sujetos: Los que han tenido dudas en algún caso frente a los que no lo han tenido en ninguno:
La característica 1 estaba presente en todos. La 4,5, 7, 8, 10 en ninguno. Por tanto solo la 2,3,6,9,11 podrían ser útiles para diferenciar a los actores/características que podrían tener algún valor para saber qué características son más fáciles de descubrir. Tomo nota para cuando nos haga falta e este vector:
La más curiosa es la característica 9. Cuando no está es fácil de acertar con ella. Cuando sí está, la tasa de acierto baja al 66%. Algo similar pasa con la 6, pero menos marcado. En todo lo demás, esté o no la característica presente es igual de fácil acertarla.
Para el resto no influye en la tasa de acierto el que esté presente o que no.
Análisis de triajes
Análisis actor a actor
En primer lugar comenzamos separando las respuestas por actores, ya que cada uno reflejaba aspectos diferentes y deberían objetivamente tener un valor de triaje diferente.
Realmente estos actores tienen una puntuación objetiva que los sujetos deberían haber asignado:
# A tibble: 7 × 2
ACTOR TRIAJEOBJ
<chr> <dbl>
1 ARTURO 2
2 MARIO 2
3 LAURA 1
4 JAVIER 1
5 MARIA 3
6 NICOLAS 1
7 MATEO 3
Vamos a evaluar los triajes que han hecho los individuos con respecto a los triajes objetivos y vamos a anotar las diferencias (si se han pasado, quedado cortos o acertado):
Me gustaría ver comparar los aciertos según GRUPO según el ACTOR, que sabemos que está muy relacionado algunos de ellos con ciertas características que estaban presentes en los vídeos:
Los aciertos corresponden al cero (VIOLETA), que desde luego no es tan abundante como esperaba. En ocasiones aciertan incluso peor que eligiendo al azar (vease el caso de ARTURO)
Los valores positivos indican que el Sujeto se ha pasado en el valor del triaje con respecto al valor objetivo, y los negativos representan lo contrario. Solo en el caso de MARIO la gente se equivoca de forma simétrica. En los demás (salvo en el de NICOLAS) todo el mundo que se equivoca, o bien se pasa o bien se queda corto. En el caso de NICOLAS, la gente tiende a quedarse corta aunque un 6% se pasa. Estaría bien pensar sobre esto.
Me resultan interesantes algunas diferencias claras a la hora de evaluar a los actores en un grupo y otro. Vamos a centrarnos en las comparaciones, actor a actor:
listaPorACTOR= dfCompTriaje %>%split(.$ACTOR) listaGraficos = listaPorACTOR%>%map2(names(listaPorACTOR), ~ {.x %>% ggstatsplot::ggbarstats(y=GRUPO, x=DIFERENCIA,bf.message =FALSE)+ggtitle(.y)})#Crear con la lista de Gráficos uno solo de dos columnaslistaGraficos
$ARTURO
$JAVIER
$LAURA
$MARIA
$MARIO
$MATEO
$NICOLAS
El único caso que queda donde es GRUPO 1 es escandalosamente peor que el GRUPO 0 es el de ARTURO. En los demás o bien el GRUPO 1 es significativamente mejor que el GRUPO 0, o bien no hay diferencia significativas entre ambos.
listaPorACTOR= dfCompTriaje %>%filter(DIFERENCIA!=0) %>%split(.$ACTOR) listaGraficos = listaPorACTOR%>%map2(names(listaPorACTOR), ~ {.x %>% ggstatsplot::ggbarstats(y=GRUPO, x=DIFERENCIA,bf.message =FALSE)+ggtitle(.y)})#Crear con la lista de Gráficos uno solo de dos columnaslistaGraficos
$ARTURO
$JAVIER
$LAURA
$MARIA
$MARIO
$MATEO
$NICOLAS
Comparaciones de resultados totales en triajes
Ahora solo nos interesa saber qué resultado total hemos tenido en los triajes. Se puede evaluar de dos formas en principio: Mirando aciertos/errores o bien viendo el error medio (en valor absoluto) cometido donde error 0 es acierto.
The Welch Two Sample t-test, which tested the difference in “Percentage of correct answers” by group (mean in intervention group = 52.56, mean in control group = 43.32), suggests that there is a positive effect from filling out the questionnaire. This effect is statistically significant and of medium magnitude (difference = 9.24, 95% CI [1.34, 17.15], t(97.22) = 2.32, p = .022; Hedge’s g = 0.44, 95% CI [0.06, 0.81]).
The Welch Two Sample t-test, which tested the difference in “mean absolute error in triaje” by group (mean in intervention group = 0.54, mean in control group = 0.68), suggests that there is a positive effect from filling out the questionnaire. This effect is statistically significant and of medium magnitude (difference = 0.18, 95% CI [0.04, 0.24], t(95.17) = 2.74, p = .007; Hedge’s g = 0.52, 95% CI [0.14, 0.89]).
Análisis multivariante
El análisis lo podemos hacer usando modelos lineales o modelos multinivel. Vamos a empezar por el lineal que siempre es más fácil de explicar.
modeloPctAciertos=lm(PCTACIERTOS ~ GRUPO+SEXO+EDAD+TMMS_ATENCION+TMMS_CLARIDAD+TMMS_REGULACION+RESILIENCIA+PSICOLOGIA+TOTALCOMP, data = dfLineal)summary(modeloPctAciertos)
Call:
lm(formula = PCTACIERTOS ~ GRUPO + SEXO + EDAD + TMMS_ATENCION +
TMMS_CLARIDAD + TMMS_REGULACION + RESILIENCIA + PSICOLOGIA +
TOTALCOMP, data = dfLineal)
Residuals:
Min 1Q Median 3Q Max
-49.912 -14.506 0.163 15.971 40.437
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 44.82081 18.99726 2.359 0.02021 *
GRUPO 10.13937 4.07237 2.490 0.01440 *
SEXO -2.35957 4.38358 -0.538 0.59156
EDAD 0.07838 0.17849 0.439 0.66148
TMMS_ATENCION 0.83306 0.59696 1.395 0.16590
TMMS_CLARIDAD 0.41879 0.71000 0.590 0.55660
TMMS_REGULACION 0.77740 0.65511 1.187 0.23811
RESILIENCIA -1.25687 0.43514 -2.888 0.00473 **
PSICOLOGIA 5.38721 4.20513 1.281 0.20306
TOTALCOMP 0.16615 0.35436 0.469 0.64016
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 20.45 on 102 degrees of freedom
(2 observations deleted due to missingness)
Multiple R-squared: 0.1494, Adjusted R-squared: 0.07436
F-statistic: 1.991 on 9 and 102 DF, p-value: 0.04778
Solo lo utilizaría para decir que se mantiene el efecto de la variable principal (incluso mejora) tras ajustar por el resto de variables. Dado que R2 no es una maravilla no me metería en intentar crear una teoría genial para ver como influye el resto de cosas en el triaje. El grupo va sobrado aúnque post-hoc podríamos estudiar un modelo donde solo esté el GRUPO y la resiliencia que es la única variable aparentemente interesante.
Los residuos no tienen mala pinta, así que por mí el modelo no está mal utilizado:
residuals(modeloPctAciertos) %>%hist()
residuals(modeloPctAciertos) %>%qqnorm()
y ahora con pruebas de normalidad de Shapiro-Wilk:
shapiro.test(residuals(modeloPctAciertos))
Shapiro-Wilk normality test
data: residuals(modeloPctAciertos)
W = 0.98579, p-value = 0.285
Vamos a hacer lo mismo para la variable ERRORMEDIO:
modeloErrorMedio=lm(ERRORMEDIO ~ GRUPO+SEXO+EDAD+TMMS_ATENCION+TMMS_CLARIDAD+TMMS_REGULACION+RESILIENCIA+PSICOLOGIA+TOTALCOMP, data = dfLineal)summary(modeloErrorMedio)
Call:
lm(formula = ERRORMEDIO ~ GRUPO + SEXO + EDAD + TMMS_ATENCION +
TMMS_CLARIDAD + TMMS_REGULACION + RESILIENCIA + PSICOLOGIA +
TOTALCOMP, data = dfLineal)
Residuals:
Min 1Q Median 3Q Max
-0.51404 -0.19494 -0.01696 0.19174 0.55809
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.7457483 0.2368059 3.149 0.00215 **
GRUPO -0.1532245 0.0507631 -3.018 0.00321 **
SEXO 0.0163033 0.0546425 0.298 0.76603
EDAD -0.0003527 0.0022249 -0.159 0.87435
TMMS_ATENCION -0.0088059 0.0074413 -1.183 0.23941
TMMS_CLARIDAD -0.0048054 0.0088504 -0.543 0.58835
TMMS_REGULACION -0.0101145 0.0081662 -1.239 0.21834
RESILIENCIA 0.0158903 0.0054241 2.930 0.00419 **
PSICOLOGIA -0.0805073 0.0524180 -1.536 0.12767
TOTALCOMP -0.0050795 0.0044172 -1.150 0.25286
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.2549 on 102 degrees of freedom
(2 observations deleted due to missingness)
Multiple R-squared: 0.1779, Adjusted R-squared: 0.1054
F-statistic: 2.453 on 9 and 102 DF, p-value: 0.01427
De nuevo, sale aún mejor que el modelo univariante. Sobre los residuos, una maravilla:
residuals(modeloErrorMedio) %>%hist()
residuals(modeloErrorMedio) %>%qqnorm()
Y ahora con pruebas de normalidad de Shapiro-Wilk
shapiro.test(residuals(modeloErrorMedio))
Shapiro-Wilk normality test
data: residuals(modeloErrorMedio)
W = 0.98886, p-value = 0.49
Modelo multinivel
EN este usamos otra base de datos. Consideramos a los actores como miembros de un conjunto infinito de actores random que van a ser evaluados por los sujetos. Con ellos podemos acertar o no en el triaje.
Ahora hacemos el modelo de regresión logístico multinivel usando el ACTOR como un factor RANDOM:
modeloMultinivel=glmer(ACIERTO ~ GRUPO+SEXO+EDAD+TMMS_ATENCION+TMMS_CLARIDAD+TMMS_REGULACION+RESILIENCIA+PSICOLOGIA+TOTALCOMP+(1|ACTOR), data = dfMultinivel, family =binomial(link ="logit"))
Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge with max|grad| = 0.00243674 (tol = 0.002, component 1)
summary(modeloMultinivel)
Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: binomial ( logit )
Formula: ACIERTO ~ GRUPO + SEXO + EDAD + TMMS_ATENCION + TMMS_CLARIDAD +
TMMS_REGULACION + RESILIENCIA + PSICOLOGIA + TOTALCOMP + (1 | ACTOR)
Data: dfMultinivel
AIC BIC logLik deviance df.resid
1047.5 1098.5 -512.8 1025.5 750
Scaled residuals:
Min 1Q Median 3Q Max
-1.9331 -0.9344 -0.6105 0.9677 1.6244
Random effects:
Groups Name Variance Std.Dev.
ACTOR (Intercept) 0.09735 0.312
Number of obs: 761, groups: ACTOR, 7
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.312940 0.745462 -0.420 0.67464
GRUPO 0.462664 0.156860 2.950 0.00318 **
SEXO -0.095900 0.170440 -0.563 0.57366
EDAD 0.003549 0.006871 0.516 0.60553
TMMS_ATENCION 0.031003 0.023089 1.343 0.17935
TMMS_CLARIDAD 0.016761 0.027397 0.612 0.54068
TMMS_REGULACION 0.029479 0.025259 1.167 0.24319
RESILIENCIA -0.048002 0.017221 -2.787 0.00531 **
PSICOLOGIA 0.155804 0.162226 0.960 0.33685
TOTALCOMP 0.008950 0.013706 0.653 0.51376
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) GRUPO SEXO EDAD TMMS_A TMMS_C TMMS_R RESILI PSICOL
GRUPO -0.020
SEXO -0.418 -0.096
EDAD -0.610 -0.236 0.258
TMMS_ATENCI -0.209 0.085 -0.314 0.244
TMMS_CLARID 0.016 -0.086 0.037 -0.063 -0.289
TMMS_REGULA -0.131 0.005 0.119 0.171 -0.008 -0.147
RESILIENCIA -0.210 0.052 0.073 -0.087 -0.098 -0.146 -0.459
PSICOLOGIA 0.117 0.090 -0.005 -0.113 0.093 -0.095 0.053 -0.104
TOTALCOMP -0.398 0.000 -0.030 0.062 -0.002 -0.275 -0.068 -0.212 -0.144
optimizer (Nelder_Mead) convergence code: 0 (OK)
Model failed to converge with max|grad| = 0.00243674 (tol = 0.002, component 1)
Los resultados son del mismo tipo, pero ahora la significación es aún mejor. Así que lo dejaría en el modelo lineal y si el revisor pide otra cosa, pues ya sabes que se lo puedes dar y saldrá aún más sólido.
La tabla ahora vendría como estimaciones de ODDS RATIOS para acertar donde solo es de interés la linea de GRUPO y el resto son solo variables de control que no necesitáis ni mostrar.
After adjusting by SEXO,…PSICOLOGIA, the odds of correctly triaging a patient in the intervention group are 1.63 times higher than in the control group (OR = 1.63, 95% CI [1.2, 2.22], p = .001).
La cuestión de evaluar cómo afecta a un profesional a cuestionar sobre su diagnóstico
La aproximación simple estaría en las pruebas chi-cuadrdo que vimos arriba. Así que si queréis hacerlo de forma más sofisticada, podríamos hacer un modelo multinivel con la variable NOTRIAJE como variable dependiente. El ACTOR seguiría siendo el efecto random y el resto variables de control con efecto fijo:
# A tibble: 2 × 5
term estimate p.value conf.low conf.high
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 0.00552 2.26e-13 0.00138 0.0222
2 GRUPO 18.4 6.20e- 5 4.43 76.6
Ahora ajustamos el modelo por el resto de variables de control:
modeloMultinivel3=glmer(NOTRIAJE ~ GRUPO+SEXO+EDAD+TMMS_ATENCION+TMMS_CLARIDAD+TMMS_REGULACION+RESILIENCIA+PSICOLOGIA+TOTALCOMP+(1|ACTOR), data = dfDuda, family =binomial(link ="logit"))
boundary (singular) fit: see help('isSingular')
summary(modeloMultinivel3)
Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: binomial ( logit )
Formula: NOTRIAJE ~ GRUPO + SEXO + EDAD + TMMS_ATENCION + TMMS_CLARIDAD +
TMMS_REGULACION + RESILIENCIA + PSICOLOGIA + TOTALCOMP + (1 | ACTOR)
Data: dfDuda
AIC BIC logLik deviance df.resid
300.5 352.1 -139.2 278.5 794
Scaled residuals:
Min 1Q Median 3Q Max
-0.6009 -0.2669 -0.1314 -0.0605 15.5069
Random effects:
Groups Name Variance Std.Dev.
ACTOR (Intercept) 6.722e-16 2.593e-08
Number of obs: 805, groups: ACTOR, 7
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -6.982207 1.924431 -3.628 0.000285 ***
GRUPO 2.797572 0.739917 3.781 0.000156 ***
SEXO 1.007049 0.435963 2.310 0.020891 *
EDAD 0.017078 0.016324 1.046 0.295478
TMMS_ATENCION 0.063177 0.051858 1.218 0.223122
TMMS_CLARIDAD -0.095585 0.063641 -1.502 0.133109
TMMS_REGULACION -0.087366 0.059498 -1.468 0.142000
RESILIENCIA 0.025824 0.038979 0.663 0.507641
PSICOLOGIA 0.981229 0.357185 2.747 0.006012 **
TOTALCOMP -0.007835 0.030132 -0.260 0.794850
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) GRUPO SEXO EDAD TMMS_A TMMS_C TMMS_R RESILI PSICOL
GRUPO -0.286
SEXO -0.528 -0.070
EDAD -0.687 -0.080 0.389
TMMS_ATENCI -0.097 0.047 -0.210 0.050
TMMS_CLARID 0.135 -0.052 -0.135 -0.193 -0.355
TMMS_REGULA -0.163 0.015 0.291 0.262 -0.236 -0.160
RESILIENCIA -0.253 0.055 0.011 0.071 0.037 -0.228 -0.416
PSICOLOGIA 0.006 0.055 -0.028 0.119 -0.017 -0.143 0.055 0.092
TOTALCOMP -0.355 -0.093 0.035 0.108 -0.004 -0.128 -0.062 -0.269 -0.287
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
EN este CASO el SEXO=1 es un factor que hace dudar y también el saber de PSICOLOGÏA.
