Ejercicio 1

bibliografía:
https://rpubs.com/gmunoz/rm-anova

library(openxlsx)
df = read.xlsx("https://docs.google.com/spreadsheets/d/e/2PACX-1vTSnSZS9JeD0AjKOwD3pNECLaQODO7T0u4E8xKi6_pts5QOAY5R_GtXDiqEl5Ta83Q-e8s0zk0Z7G1G/pub?output=xlsx",sheet="16_anova_repetidas")
head(df)
##   ID Values Month
## 1  1    108 mes30
## 2  2    103 mes30
## 3  3     96 mes30
## 4  4     84 mes30
## 5  5    118 mes30
## 6  6    110 mes30

convertir a factores, como son menos de 50 observaciones podemos usar shapiro

shapiro.test(df$Values[df$Month=="mes30"])
## 
##  Shapiro-Wilk normality test
## 
## data:  df$Values[df$Month == "mes30"]
## W = 0.96751, p-value = 0.883
df$ID<-factor(df$ID)
df$Month<-factor(df$Month)
str(df)
## 'data.frame':    48 obs. of  3 variables:
##  $ ID    : Factor w/ 12 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ Values: num  108 103 96 84 118 110 129 90 84 96 ...
##  $ Month : Factor w/ 4 levels "mes30","mes36",..: 1 1 1 1 1 1 1 1 1 1 ...

Necesitamos ver la normalidad por grupo

tapply( df$Values,df$Month, shapiro.test)
## $mes30
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.96751, p-value = 0.883
## 
## 
## $mes36
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.95671, p-value = 0.7361
## 
## 
## $mes42
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.91137, p-value = 0.2222
## 
## 
## $mes48
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.95958, p-value = 0.7778

Se cumple la normalidad en todos los casos

Primer método

No da esfericidad

rm.aov<-aov(Values~Month+Error(ID/Month),data=df)
summary(rm.aov)
## 
## Error: ID
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 11   6624   602.2               
## 
## Error: ID:Month
##           Df Sum Sq Mean Sq F value Pr(>F)  
## Month      3    552  184.00   3.027 0.0432 *
## Residuals 33   2006   60.79                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Segundo método

No da esfericidad

demo1.aov <- aov(Values ~ Month * Month + Error(ID), data = df)
summary(demo1.aov)
## 
## Error: ID
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 11   6624   602.2               
## 
## Error: Within
##           Df Sum Sq Mean Sq F value Pr(>F)  
## Month      3    552  184.00   3.027 0.0432 *
## Residuals 33   2006   60.79                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Tercer método

library(ez)
#Aplicamos un medidas repetidas:
ezANOVA(data=df, dv=Values, wid=ID, within=Month)
## $ANOVA
##   Effect DFn DFd        F          p p<.05        ges
## 2  Month   3  33 3.026919 0.04321863     * 0.06011762
## 
## $`Mauchly's Test for Sphericity`
##   Effect         W          p p<.05
## 2  Month 0.2426472 0.01771762     *
## 
## $`Sphericity Corrections`
##   Effect       GGe     p[GG] p[GG]<.05       HFe      p[HF] p[HF]<.05
## 2  Month 0.6095445 0.0747874           0.7248502 0.06353773

La esfericidad es menor de 0.05 entonces no cumple cumple con la homogeneidad de varianzas ni tampoco ayudan las correcciones de esfericidad ya que estas son mayores a 0.056, se debe realizar una prueba no paramétrica

Friedman

friedman.test(df$Values, df$Month, df$ID)
## 
##  Friedman rank sum test
## 
## data:  df$Values, df$Month and df$ID
## Friedman chi-squared = 5.819, df = 3, p-value = 0.1208

no hay diferencias estadísticas

pairwise.wilcox.test(df$Values, df$Month, paired = TRUE, p.adjust.method = "holm")
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with zeroes

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with zeroes
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with zeroes
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## 
##  Pairwise comparisons using Wilcoxon signed rank test with continuity correction 
## 
## data:  df$Values and df$Month 
## 
##       mes30 mes36 mes42
## mes36 0.63  -     -    
## mes42 0.25  1.00  -    
## mes48 0.25  1.00  1.00 
## 
## P value adjustment method: holm

##Ejercicio 2

Buena bibliografia:
https://www.maximaformacion.es/blog-dat/realiza-e-interpreta-paso-a-paso-un-anova-de-medidas-repetidas-en-r-software/

library(WRS2)
data("hangover")
head(hangover)
##   symptoms   group time id
## 1        0 control    1  1
## 2       32 control    1  2
## 3        9 control    1  3
## 4        0 control    1  4
## 5        2 control    1  5
## 6        0 control    1  6
str(hangover)
## 'data.frame':    120 obs. of  4 variables:
##  $ symptoms: num  0 32 9 0 2 0 41 0 0 0 ...
##  $ group   : Factor w/ 2 levels "control","alcoholic": 1 1 1 1 1 1 1 1 1 1 ...
##  $ time    : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
##  $ id      : int  1 2 3 4 5 6 7 8 9 10 ...
hangover$id<-factor(hangover$id) #transformamos id a factor

Vemos que tenemos 4 variables:

symptoms: Número de síntomas de la resaca
group: Hijo de alcohólico vs. control
time: Momento de la medición
ID: identificación del sujeto
Aquí no utilizaremos la variable “group” que diferencia entre hijos de padres alcohólicos y los que no, pero podríamos incluir este factor de grupo, y su posible interacción con el tiempo de medición, utilizando un ANOVA mixto (lo dejo para vosotros).

Una vez tenemos los datos, vamos a realizar un análisis descriptivo para tener una idea de cómo son nuestros datos y luego analizamos los supuestos de la prueba 

by(hangover$symptoms,hangover$time, length)
## hangover$time: 1
## [1] 40
## ------------------------------------------------------------ 
## hangover$time: 2
## [1] 40
## ------------------------------------------------------------ 
## hangover$time: 3
## [1] 40

Vemos que tenemos 40 casos (o sujetos) para cada una de las 3 muestras de tiempo.

by(hangover$symptoms,hangover$time, summary)
## hangover$time: 1
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   0.000   5.175   6.500  41.000 
## ------------------------------------------------------------ 
## hangover$time: 2
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   4.000   6.875  11.250  43.000 
## ------------------------------------------------------------ 
## hangover$time: 3
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   3.000   6.625  10.250  42.000

Realizamos ahora un resumen numérico para cada muestra y el gráfico de cajas (boxplot):

boxplot(hangover$symptoms~hangover$time, col=c("lightblue","pink","purple"))

Vemos que no varía demasiado el número de síntomas en las 3 mediciones, aunque el tiempo 1 parece tener menor número existe mucha dispersión en los datos. También observamos la presencia de datos atípicos (outliers; los puntos del gráfico) que podrían indicarnos que es mejor realizar pruebas robustas para el ANOVA MR.

Ahora comprobamos el supuesto de normalidad.

by(hangover$symptoms,hangover$time, shapiro.test)
## hangover$time: 1
## 
##  Shapiro-Wilk normality test
## 
## data:  dd[x, ]
## W = 0.59901, p-value = 3.033e-09
## 
## ------------------------------------------------------------ 
## hangover$time: 2
## 
##  Shapiro-Wilk normality test
## 
## data:  dd[x, ]
## W = 0.76198, p-value = 1.195e-06
## 
## ------------------------------------------------------------ 
## hangover$time: 3
## 
##  Shapiro-Wilk normality test
## 
## data:  dd[x, ]
## W = 0.75338, p-value = 8.256e-07

En los 3 tiempos la variable “symptoms” no tiene distribución normal, sería interesante desarrollar pruebas no paramétricas que no trabajan bien para datos que no son normales.

También podríamos realizar un gráfico Q-Q para evaluar este supuesto (lo dejo en vuestras manos!).

par(mfrow=c(2,2))

for (i in 1:3){
  hist(hangover[hangover$time==i,1], main=("Time"), 
       breaks = 20)
}
par(mfrow=c(1,1))

Se necesita una prueba de friedman o probar con eliminar los valores outliers,

str(hangover)
## 'data.frame':    120 obs. of  4 variables:
##  $ symptoms: num  0 32 9 0 2 0 41 0 0 0 ...
##  $ group   : Factor w/ 2 levels "control","alcoholic": 1 1 1 1 1 1 1 1 1 1 ...
##  $ time    : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
##  $ id      : Factor w/ 40 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...

si lo hacemos con la prueba parametrica ez anova, no son significativos

friedman.test(hangover$symptoms, hangover$time, hangover$id)
## 
##  Friedman rank sum test
## 
## data:  hangover$symptoms, hangover$time and hangover$id
## Friedman chi-squared = 8.887, df = 2, p-value = 0.01175

Si hay diferencias significativas en los 3 tiempos p=0.01175, es decir, el número de síntomas de resaca si difiere significativamente entre los 3 tiempos de medición.

pairwise.wilcox.test(hangover$symptoms, hangover$time, paired = TRUE, p.adjust.method = "holm")
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with zeroes
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with zeroes
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with zeroes
## 
##  Pairwise comparisons using Wilcoxon signed rank test with continuity correction 
## 
## data:  hangover$symptoms and hangover$time 
## 
##   1    2   
## 2 0.21 -   
## 3 0.21 0.98
## 
## P value adjustment method: holm

marca un error, hay que checarlo

Ejercicio Nuevo

bibliografia buena
https://rstudio-pubs-static.s3.amazonaws.com/511966_05aefff3aaca41818010b40aa05bfece.html


Queremos realizar un ANÁLISIS DE MEDIDAS REPETIDAS para los datos de rendimiento académico de seis alumnos. En cada alumno se ha medido a cinco tiempos diferentes su rendimiento, por tanto las muestras tomadas no son independientes entre si. Para poder analizar estos datos debemos considerar las muestras como relacionadas, es decir debemos realizar un ANOVA de medidas repetidas.

Veamos como son nuestros datos:

datosmedidas = read.xlsx("https://docs.google.com/spreadsheets/d/e/2PACX-1vTSnSZS9JeD0AjKOwD3pNECLaQODO7T0u4E8xKi6_pts5QOAY5R_GtXDiqEl5Ta83Q-e8s0zk0Z7G1G/pub?output=xlsx",sheet="19_rep")
head(datosmedidas)
##   individuo time rendimiento
## 1         1    1         8.5
## 2         1    2         8.2
## 3         1    3         8.9
## 4         1    4         7.7
## 5         1    5         7.4
## 6         2    1         9.8
datosmedidas$time<-factor(datosmedidas$time)
datosmedidas$individuo<-factor(datosmedidas$individuo)
library(ez)
ezANOVA(data=datosmedidas, dv=rendimiento, wid=individuo, within=time)
## $ANOVA
##   Effect DFn DFd        F            p p<.05       ges
## 2   time   4  20 12.41317 3.096979e-05     * 0.2211562
## 
## $`Mauchly's Test for Sphericity`
##   Effect           W          p p<.05
## 2   time 0.003212948 0.03413309     *
## 
## $`Sphericity Corrections`
##   Effect       GGe       p[GG] p[GG]<.05       HFe       p[HF] p[HF]<.05
## 2   time 0.3126606 0.009730034         * 0.3669864 0.006076078         *

Podemos ver en nuestros resultados que el tiempo es significativo ($ANOVA el p es menor de 0.05), es decir el rendimiento escolar cambia con el tiempo. Pero en este caso nuestros datos violan uno de los requisitos que es la esfericidad (Mauchly’s test es significativo), debemos fiarnos de la p de la ‘Sphericity corrections’ que nos confirma lo que hemos deducido al principio.

boxplot(rendimiento~time, xlab="time", ylab="Rendimiento académico", main="rendimiento alumnos con el paso del tiempo", col="lightblue", data=datosmedidas)

Claramente el tiempo hace que el rendimiento académico decaiga, ¿se aburren nuestros alumnos?

Si queremos ver el rendimiento de cada uno de nuestros alumnos:


Prueba no paramétrica

Nuestros datos deben estar en otro orden

library(openxlsx)
datosmedidasc = read.xlsx("https://docs.google.com/spreadsheets/d/e/2PACX-1vTSnSZS9JeD0AjKOwD3pNECLaQODO7T0u4E8xKi6_pts5QOAY5R_GtXDiqEl5Ta83Q-e8s0zk0Z7G1G/pub?output=xlsx",sheet="20_rep")
head(datosmedidasc)
##   individuo time1 time2 time3 time4 time5
## 1         1   8.5   8.2   8.9   7.7   7.4
## 2         2   9.8   8.9   8.9   8.8   8.1
## 3         3   9.6   9.0   9.3   7.5   7.1
## 4         4   7.5   7.8   7.8   4.5   4.6
## 5         5   5.8   5.8   5.9   2.6   1.2
## 6         6   9.9   9.8   9.6   8.6   8.7
library(jmv)
ann<-anovaRMNP(datosmedidasc,measures=vars(time1,time2,time3,time4,time5))
ann
## 
##  REPEATED MEASURES ANOVA (NON-PARAMETRIC)
## 
##  Friedman                        
##  ------------------------------- 
##    <U+03C7>²          df    p           
##  ------------------------------- 
##    19.24786     4    0.0007025   
##  -------------------------------
ann$table
## 
##  Friedman                        
##  ------------------------------- 
##    <U+03C7>²          df    p           
##  ------------------------------- 
##    19.24786     4    0.0007025   
##  -------------------------------

La conclusión es la misma, el tiempo influye en el rendimiento académico de los alumnos.

Ejercicio Nuevo

bibliografia: https://www.youtube.com/watch?v=OVrdl-j_pk8&t=1124s

base = read.xlsx("https://docs.google.com/spreadsheets/d/e/2PACX-1vTSnSZS9JeD0AjKOwD3pNECLaQODO7T0u4E8xKi6_pts5QOAY5R_GtXDiqEl5Ta83Q-e8s0zk0Z7G1G/pub?output=xlsx",sheet="21_rep")
head(base)
##   id   tiempo puntuacion
## 1  1 tiempo_1         61
## 2  2 tiempo_1         64
## 3  3 tiempo_1         62
## 4  4 tiempo_1         60
## 5  5 tiempo_1         65
## 6  6 tiempo_1         64
str(base)
## 'data.frame':    30 obs. of  3 variables:
##  $ id        : num  1 2 3 4 5 6 7 8 9 10 ...
##  $ tiempo    : chr  "tiempo_1" "tiempo_1" "tiempo_1" "tiempo_1" ...
##  $ puntuacion: num  61 64 62 60 65 64 63 65 60 61 ...
base$tiempo<-factor(base$tiempo)
base$id<-factor(base$id)
boxplot(base$puntuacion~base$tiempo, 
        col=c("purple","orange","pink")
        )

#tapply( df$puntuacion,df$tiempo, length)

menos de 30 observaciones

#tapply( df$puntuacion,df$tiempo, shapiro.test)

se distribuyen de manera normal

Hacer modelo

library(rstatix)

modelo1<-anova_test(
  data=base,
  dv=puntuacion,
  wid=id,
  within=tiempo)
modelo1$ANOVA
##   Effect DFn DFd       F        p p<.05   ges
## 1 tiempo   2  18 173.162 1.75e-12     * 0.914

si hay diferencias a menos en un grupo

modelo1$`Mauchly's Test for Sphericity`
##   Effect     W    p p<.05
## 1 tiempo 0.823 0.46

cumplen con el supuesto de esfericidad (varianzas iguales)

modelo1$`Sphericity Corrections`
##   Effect  GGe    DF[GG]    p[GG] p[GG]<.05   HFe      DF[HF]    p[HF] p[HF]<.05
## 1 tiempo 0.85 1.7, 15.3 6.74e-11         * 1.027 2.05, 18.49 1.75e-12         *

Poner las correcciones por si no se cumple la esfericidad 

#agregamos el pes (eta cuadrada parcial)
modelo1<-anova_test(
  data=base,
  dv=puntuacion,
  wid=id,
  within=tiempo,
  effect.size = "pes"
  )
modelo1$ANOVA
##   Effect DFn DFd       F        p p<.05   pes
## 1 tiempo   2  18 173.162 1.75e-12     * 0.951

el 95% de los cambios se explica por los tratamientos

comparaciones

tiempo1Valores<-base[1:10,3]
tiempo2Valores<-base[11:20,3]
tiempo3Valores<-base[21:30,3]
t1<-t.test(tiempo1Valores,tiempo2Valores,paired=T)
t2<-t.test(tiempo1Valores,tiempo3Valores,paired=T)
t3<-t.test(tiempo2Valores,tiempo3Valores,paired=T)
#Report the pvalues of each test in one row.
pvalues<-c(t1$p.value,t2$p.value,t3$p.value)
#Adjust the p-values using the bonferroni correction for 3 tests (rounded to 4 decimals).
round(p.adjust(pvalues,'bonferroni',3),10) #bonferroni o holm
## [1] 0.0000000854 0.0318413207 0.0000016539

Tres ejercicios mas para checar

# https://personal.us.es/vararey/adatos2/materiales/anovarepe.pdf
# https://rpubs.com/Joaquin_AR/219732
# https://www.youtube.com/watch?v=UejePH5dAkQ