── 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)
Warning: package 'ggstatsplot' was built under R version 4.4.1
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)
Warning: package 'lme4' was built under R version 4.4.1
Loading required package: Matrix
Attaching package: 'Matrix'
The following objects are masked from 'package:tidyr':
expand, pack, unpack
library(effects)
Warning: package 'effects' was built under R version 4.4.1
Loading required package: carData
lattice theme set by effectsTheme()
See ?effectsTheme for details.
library(ggeffects)
Warning: package 'ggeffects' was built under R version 4.4.1
library(emmeans)
Warning: package 'emmeans' was built under R version 4.4.1
Welcome to emmeans.
Caution: You lose important information if you filter this package's results.
See '? untidy'
library(gtsummary)
Warning: package 'gtsummary' was built under R version 4.4.1
dfRaw <-read_sav("TRIAGE_REDUCIDO.sav") %>%# convertir la columna grupo en factor, con los niveles 0 y 1 convertidos en "Without checklist", "With checklist"mutate(Group =factor(as.factor(GRUPO), levels =c(0, 1), labels =c("Without checklist", "With checklist"))) %>%mutate(SEXO=factor(as.factor(SEXO), levels =c(1, 2), labels =c("Hombre", "Mujer"))) %>%mutate(PSICOLOGIA=factor(as.factor(PSICOLOGIA), levels =c(0, 1), labels =c("Otro", "Psicólogo")))
Hay variables que se usarán posteriormente y las aparto ahora:
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:
Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.
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:
#unir gráficos p1 y p2 en uno sologridExtra::grid.arrange(p1,p2, nrow=1)
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:
Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.
p3
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.
Vamos a desglosar esos patrones según GRUPO:
p3<- dfCompTriaje %>%filter(GRUPO==0) %>%ggbarstats(y=Actor, x=`Triage results`, bf.message =FALSE, title="Cases without checklist")+theme(axis.text.x =element_text(angle =90, hjust =1, vjust =0.5))+xlab("")+# usar paleta ordinal para el colorscale_fill_manual(values =c("darkblue", "blue", "gray", "orange", "darkorange"))
Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.
p3
p4<- dfCompTriaje %>%filter(GRUPO==1) %>%ggbarstats(y=Actor, x=`Triage results`, bf.message =FALSE, title="Cases with checklist")+theme(axis.text.x =element_text(angle =90, hjust =1, vjust =0.5))+xlab("")+# usar paleta ordinal para el colorscale_fill_manual(values =c("darkblue", "blue", "gray", "orange", "darkorange"))
Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.
p4
Los dos juntos:
#unir gráficos p3 y p4 en uno sologridExtra::grid.arrange(p3,p4, nrow=2)
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.
dfTriajeTotal %>%ggbetweenstats(x=Group, y=ERRORMEDIO, bf.message =FALSE)+ylab("Error medio cometido")+xlab("")
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 ~ Group+SEXO+EDAD+TMMS_ATENCION+TMMS_CLARIDAD+TMMS_REGULACION+RESILIENCIA+PSICOLOGIA+#TOTALCOMP+ (1|ACTOR), data = dfMultinivel, family =binomial(link ="logit"))summary(modeloMultinivel)
Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: binomial ( logit )
Formula: ACIERTO ~ Group + SEXO + EDAD + TMMS_ATENCION + TMMS_CLARIDAD +
TMMS_REGULACION + RESILIENCIA + PSICOLOGIA + (1 | ACTOR)
Data: dfMultinivel
AIC BIC logLik deviance df.resid
1046.0 1092.3 -513.0 1026.0 751
Scaled residuals:
Min 1Q Median 3Q Max
-1.9047 -0.9356 -0.6215 0.9677 1.5628
Random effects:
Groups Name Variance Std.Dev.
ACTOR (Intercept) 0.09721 0.3118
Number of obs: 761, groups: ACTOR, 7
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.211819 0.622685 -0.340 0.73373
GroupWith checklist 0.462874 0.156869 2.951 0.00317 **
SEXOMujer -0.092525 0.170230 -0.544 0.58676
EDAD 0.003273 0.006861 0.477 0.63332
TMMS_ATENCION 0.031009 0.023070 1.344 0.17890
TMMS_CLARIDAD 0.021772 0.026276 0.829 0.40733
TMMS_REGULACION 0.030604 0.025183 1.215 0.22426
RESILIENCIA -0.045663 0.016807 -2.717 0.00659 **
PSICOLOGIAPsicólogo 0.170821 0.160448 1.065 0.28703
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Warning in abbreviate(rn, minlength = 11): abbreviate used with non-ASCII chars
Correlation of Fixed Effects:
Warning in abbreviate(rn, minlength = 6): abbreviate used with non-ASCII chars
Y si quieres que las tablas tengan el número de aciertos para cada nivel de una variable y dar más detalle como los q-valores (p-valores corregidos para comparaciones múltiples) y factor de inflación de la varianza (VIF) por si queréis comentar que no había mucha colinealidad entre las variables explicativas. Más detalles que esto no se suele dar. Podéis borrar las columnas que no os gusten:
OR = Odds Ratio, CI = Confidence Interval, VIF = Variance Inflation Factor
3
False discovery rate correction for multiple testing
O así:
sjPlot::tab_model(modeloMultinivel)
ACIERTO
Predictors
Odds Ratios
CI
p
(Intercept)
0.81
0.24 – 2.74
0.734
Group: With checklist
1.59
1.17 – 2.16
0.003
SEXO: Mujer
0.91
0.65 – 1.27
0.587
EDAD
1.00
0.99 – 1.02
0.633
FACTOR ATENCIÓN INT. EMOCIONAL
1.03
0.99 – 1.08
0.179
FACTOR CLARIDAD INT. EMOCIONAL
1.02
0.97 – 1.08
0.407
FACTOR REGULACIÓN INT. EMOCIONAL
1.03
0.98 – 1.08
0.224
RESILIENCIA
0.96
0.92 – 0.99
0.007
PSICOLOGIA: Psicólogo
1.19
0.87 – 1.62
0.287
Random Effects
σ2
3.29
τ00ACTOR
0.10
ICC
0.03
N ACTOR
7
Observations
761
Marginal R2 / Conditional R2
0.032 / 0.060
Una forma de mostrar gráficamente qué significa la parte interesante del modelo es centrarnos en el efecto del grupo en la probabilidad de acierto. Para ello, podemos usar la función ggpredict del paquete ggeffects para generar predicciones.
# Generar predicciones usando ggpredictpredicciones <-ggpredict(modeloMultinivel, terms ="Group")# Graficar las prediccionesplot(predicciones) +theme_minimal() +labs(title ="Efecto del grupo en la probabilidad de acierto",x ="Group",y ="Probabilidad de acierto")+coord_cartesian(ylim=c(0.3,0.65))
Group prob SE df asymp.LCL asymp.UCL null z.ratio p.value
Without checklist 0.431 0.0402 Inf 0.354 0.510 0.5 -1.704 0.0883
With checklist 0.546 0.0393 Inf 0.468 0.621 0.5 1.157 0.2473
Results are averaged over the levels of: SEXO, PSICOLOGIA
Confidence level used: 0.95
Intervals are back-transformed from the logit scale
Tests are performed on the logit scale
Esos “p” de ahí tratan sobre si la probabilidad de acertar se aleja del 50%, cosa que no es relevante, pero lo pongo por si os interesa tener el intervalo de confianza de la probabilidad de acertar en forma numérica en lugar de gráfica.
Y de nuevo la OR en forma de solo valores numéricos para la variable interesante. Lo pongo así por tener un doble chequeo de que lo haga como lo haga da lo mismo.
contrast odds.ratio SE df asymp.LCL asymp.UCL
With checklist / Without checklist 1.59 0.249 Inf 1.17 2.16
null z.ratio p.value
1 2.951 0.0032
Results are averaged over the levels of: SEXO, PSICOLOGIA
Confidence level used: 0.95
Intervals are back-transformed from the log odds ratio scale
Tests are performed on the log odds ratio scale
