Primero vamos a cargar los datos del fichero spss:
library(tidyverse)
Warning: package 'ggplot2' was built under R version 4.3.2
Warning: package 'dplyr' was built under R version 4.3.2
Warning: package 'stringr' was built under R version 4.3.2
── 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.4.4 ✔ 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(haven)
Warning: package 'haven' was built under R version 4.3.2
library(ggstatsplot)
Warning: package 'ggstatsplot' was built under R version 4.3.2
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) No calcula p, asi que la cambio por otra:library(lmerTest)
Warning: package 'lmerTest' was built under R version 4.3.2
Loading required package: lme4
Warning: package 'lme4' was built under R version 4.3.2
Loading required package: Matrix
Warning: package 'Matrix' was built under R version 4.3.2
Attaching package: 'Matrix'
The following objects are masked from 'package:tidyr':
expand, pack, unpack
Attaching package: 'lmerTest'
The following object is masked from 'package:lme4':
lmer
The following object is masked from 'package:stats':
step
library(patchwork)library(gtsummary)
Warning: package 'gtsummary' was built under R version 4.3.2
library(afex)
Warning: package 'afex' was built under R version 4.3.2
************
Welcome to afex. For support visit: http://afex.singmann.science/
- Functions for ANOVAs: aov_car(), aov_ez(), and aov_4()
- Methods for calculating p-values with mixed(): 'S', 'KR', 'LRT', and 'PB'
- 'afex_aov' and 'mixed' objects can be passed to emmeans() for follow-up tests
- Get and set global package options with: afex_options()
- Set sum-to-zero contrasts globally: set_sum_contrasts()
- For example analyses see: browseVignettes("afex")
************
Attaching package: 'afex'
The following object is masked from 'package:lme4':
lmer
Nos gustaría estudiar si la variable LongBalanceo presenta diferencias entre Ojos CE y OE ajustando por inmovilización. Para ello, vamos a realizar un modelo lineal mixto con la función lmer de la librería lme4. Nos sale que sí ocurre que hay diferencias.
Queremos también saber si el salto que se produce en Inmov=No entre ojos cerrados y abiertos es diferente al de cuando hay inmovilización dominante o no dominante. Transformamos la base de datos y exploramos las primers líneas:
En forma de gráfico, el modelo sin considerar interacciones con los ojos es este (es como tener la media de ojos aboertos y cerrados y estudiar solo el efecto de la inmovilización).
p0<-ggwithinstats(dg, subject.id="DNI", x = Inmov, y = Valor, pairwise.comparisons =TRUE, pairwise.display ="significant", pairwise.grouping.line.type ="segment",bf.message=FALSE,point.args =list(position =position_jitter(width =0.1, height =0)) ) +coord_cartesian(ylim=c(0,150))+labs(title="Ignorando ojos",y="Longitud de Balanceo")p0
Para más detalle vamos a estudiarlo para experimentos con ojos abiertos y cerrados por separado que es más o menos lo que hace el modelo con interacciones.
p1<-ggwithinstats(dg %>%filter(Ojos=="CE"),subject.id="DNI", x = Inmov, y = Valor, pairwise.comparisons =TRUE, pairwise.display ="significant", pairwise.grouping.line.type ="segment",bf.message=FALSE,point.args =list(position =position_jitter(width =0.1, height =0)) ) +coord_cartesian(ylim=c(0,150))+labs(title="Closed eyes",y="Longitud de Balanceo")p1
p2<-ggwithinstats(dg %>%filter(Ojos=="OE"), x = Inmov, y = Valor, pairwise.comparisons =TRUE, pairwise.display ="significant", pairwise.grouping.line.type ="segment",bf.message=FALSE, point.args =list(position =position_jitter(width =0.1, height =0))) +coord_cartesian(ylim=c(0,150))+labs(title="Open eyes",y="Longitud de Balanceo")p2
Ahora me gustaría ensamblar los gráficos p1 y p2 en uno solo.
patchwork::wrap_plots(p1,p2,ncol=2)
Queda claro que hay un efecto de la inmovilización esten como estén los ojos. No hay diferencia entre que se haga en el miembro dominante o no dominante cn los ojos abiertos, aunque sí que hay una pequeña con los ojos cerrados. No sabría si prestarle mucha atención ya que las interacciones aunque salieron cerca de la significación, no lo fueron. Puedes elegir qué escribir.
Las comparaciones post-hoc están ajustadas para comparaciones múltiples por el método de Holm.
Todo lo anterior es dependiente de que se den las condiciones de validez (esfericidad y normalidad). Veamos si hay algo en contra de ellas:
pruebaAnova <-aov_ez(dg, id ="DNI", dv ="Valor",within =c("Ojos", "Inmov"), detailed =TRUE)pruebaAnova %>%summary()
Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
Sum Sq num Df Error SS den Df F value Pr(>F)
(Intercept) 617142 1 84042 40 293.7314 < 2.2e-16 ***
Ojos 5559 1 7171 40 31.0071 1.909e-06 ***
Inmov 5578 2 11168 80 19.9771 9.183e-08 ***
Ojos:Inmov 493 2 8300 80 2.3739 0.09965 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Mauchly Tests for Sphericity
Test statistic p-value
Inmov 0.82100 0.021365
Ojos:Inmov 0.89465 0.114098
Greenhouse-Geisser and Huynh-Feldt Corrections
for Departure from Sphericity
GG eps Pr(>F[GG])
Inmov 0.84818 6.609e-07 ***
Ojos:Inmov 0.90469 0.1055
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
HF eps Pr(>F[HF])
Inmov 0.8817775 4.267377e-07
Ojos:Inmov 0.9450613 1.029672e-01
Como la prueba de Mauchly es significativa para Inmov (en ojos no hay que hacerlo ya que tiene dos niveles y con esas no es necesario) se puede hacer la corrección de Greenhouse-Geisser. En este caso no cambia nada, así que el asunto de la esfericidad no nos preocupa en las interpretaciones.
residuals(pruebaAnova) %>%shapiro.test()
Shapiro-Wilk normality test
data: .
W = 0.98316, p-value = 0.005202
En cuanto a la normalidad hemos tenido evidencia en contra. Sin embargo no es algo que me preocupe demasiado ya que el tamaño muestral es razonablemente grande y la prueba de normalidad es muy sensible a ello. Además, el modelo mixto es algo robusto a la violación de la normalidad especialmente en casos como este donde no veo observaciones anómalas. Además, mira esto:
residuals(pruebaAnova) %>%hist()
residuals(pruebaAnova) %>%scale() %>%qqnorm()
¿Qué problema hay? En peores plazas hemos toreao.
Objetivo 2
Estudiar si la Superficie total cuando No hay Inmobilización es es más baja que cuando sí la hay, sea dominante o no dominante.
Pues nada… Aquí solo hay un efecto de los ojos. La inmovilización no juega gran cosa. Por curiosidad vamos a visualizarlo, pero no se necesita. Ademas en el siguiente gráfico se ve que los datos no tienen distribución normal, así que paso a mostrarlo con pruebas no paramétricas.
