Medidas Repetidas_1,2 y 3 vias
Palma Oliva Anderson
2023-06-02
#install.packages("datarium")
#data("selfesteem", package = "datarium")
#datos = selfesteem
#head(datos, 3)#convertir formato ancho en formato largo
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(tidyr)
# install.packages("datarium")
data("selfesteem", package = "datarium")
datos = selfesteem
datos = datos %>%
gather(key = "tiempo",
value = "rto",
t1, t2, t3) %>%
mutate_at(vars(id, tiempo), as.factor) #aqui el tiempo y el id se concirtieron a factor para que R no los tome como numeros
head(datos)## # A tibble: 6 × 3
## id tiempo rto
## <fct> <fct> <dbl>
## 1 1 t1 4.01
## 2 2 t1 2.56
## 3 3 t1 3.24
## 4 4 t1 3.42
## 5 5 t1 2.87
## 6 6 t1 2.05#View(datos)
#gather es una funcion que agrupa datos, en este caso pone todos los rendiemientos en una sola columna***Convertir todos los datos de formato ancho a formato largo Resumen estadistico de los datos
****Resumen estadistico por tabulacion
#library(rst)
#resumen estadistico de los datos
datos %>%
group_by(tiempo) %>%
summarise(media = mean(rto),
desv = sd(rto),
n = n(),
cv = 100*desv/media)## # A tibble: 3 × 5
## tiempo media desv n cv
## <fct> <dbl> <dbl> <int> <dbl>
## 1 t1 3.14 0.552 10 17.6
## 2 t2 4.93 0.863 10 17.5
## 3 t3 7.64 1.14 10 15.0#Nota: coeficientes de vatiacion menors al 20% ususlamente se asocian a datos que tienen comportamiento normal #Cuando son mayores generalmente representan problemas, como heterogeneidad, datos extraños, datos mal tomas, generalmente requiere imputacion o causa datos atipicos
*****Resumen estadistico por vusualizacion
boxplot(datos$rto ~ datos$tiempo)
Deteccion de outliers
# (deteccion de datos atipicos)
library(rstatix)##
## Attaching package: 'rstatix'## The following object is masked from 'package:stats':
##
## filter#atipicos
datos %>%
group_by(tiempo) %>%
identify_outliers(rto)## # A tibble: 2 × 5
## tiempo id rto is.outlier is.extreme
## <fct> <fct> <dbl> <lgl> <lgl>
## 1 t1 6 2.05 TRUE FALSE
## 2 t2 2 6.91 TRUE FALSE#Existen 2 datos que corren el peligro de ser atipicos (TRUE), t1 id6 y t2 id2, sin embargo si se observa el extremo este es FALSE, ademas el coeficiente de variacion fue menor al 20 % por lo que se descarta que sean datos atipicos
Normalidad
datos %>%
group_by(tiempo) %>%
shapiro_test(rto)## # A tibble: 3 × 4
## tiempo variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 t1 rto 0.967 0.859
## 2 t2 rto 0.876 0.117
## 3 t3 rto 0.923 0.380#Aqui se observa la columna p, este p-valor debe ser mayor a 5%, en este caso tanto t1, como t2 y t3 son mayores al 5%, es decir que los datos en cada tiempo son normales ##En caso de ser menor a 5% se pone una nota que diga los datos de este tiempo parecen no tener comportamiento normal, esto puede ser causada por un datos atipico que reuqiera la eliminacion del dato o la imputacion del dato
Supuestos de Esfericidad
#supuestos de esfericidad (igualdad de varianzas)
res.aov <- anova_test(data = datos,
dv = rto,
wid = id,
within = tiempo)
res.aov$`Mauchly's Test for Sphericity`## Effect W p p<.05
## 1 tiempo 0.551 0.092get_anova_table(res.aov)## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 tiempo 2 18 55.469 2.01e-08 * 0.829#p-valor=0.092 (9.2%) como es mayor al 5% puedo asumir varianzas iguales #(rconsole), p-valor=2.01e-08, como es menor al 5%, rechaza la hipotesis nula ##Como se rechaza la hipotesis nula se concluye que el aceite que se produce en los 3 tiempos no es el mismo (el corte en el tiempo si tiene efecto, en el corte 3 es el mayor rendimineto de aceite)
# se usa el p-valor ajustado (p.adj)
datos %>%
pairwise_t_test(
rto ~ tiempo,
paired = TRUE,
p.adjust.method = "bonferroni")## # A tibble: 3 × 10
## .y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 rto t1 t2 10 10 -4.97 9 0.000772 2e-3 **
## 2 rto t1 t3 10 10 -13.2 9 0.000000334 1e-6 ****
## 3 rto t2 t3 10 10 -4.87 9 0.000886 3e-3 **Se observa p.adj, del corte 1 al corte 3 el rendimineto baja
**********************************************************************+
Medidas repetidas de dos vias
data("selfesteem2", package = "datarium")
datos2 = selfesteem2
datos2$treatment = gl(2,12,24,
c('con fert',
'sin fert'))
head(datos2)## # A tibble: 6 × 5
## id treatment t1 t2 t3
## <fct> <fct> <dbl> <dbl> <dbl>
## 1 1 con fert 83 77 69
## 2 2 con fert 97 95 88
## 3 3 con fert 93 92 89
## 4 4 con fert 92 92 89
## 5 5 con fert 77 73 68
## 6 6 con fert 72 65 63datos2 = datos2 %>%
gather(key='tiempo', value = 'rto',
t1,t2,t3)
#View(datos2)
head(datos2)## # A tibble: 6 × 4
## id treatment tiempo rto
## <fct> <fct> <chr> <dbl>
## 1 1 con fert t1 83
## 2 2 con fert t1 97
## 3 3 con fert t1 93
## 4 4 con fert t1 92
## 5 5 con fert t1 77
## 6 6 con fert t1 72Resumen estadistico tabular #Cambiar los datos a formato largo
datos2 %>%
group_by(treatment, tiempo) %>%
summarise(media = mean(rto),
desv = sd(rto),
n = n(),
cv = 100*desv/media)## `summarise()` has grouped output by 'treatment'. You can override using the
## `.groups` argument.## # A tibble: 6 × 6
## # Groups: treatment [2]
## treatment tiempo media desv n cv
## <fct> <chr> <dbl> <dbl> <int> <dbl>
## 1 con fert t1 88 8.08 12 9.18
## 2 con fert t2 83.8 10.2 12 12.2
## 3 con fert t3 78.7 10.5 12 13.4
## 4 sin fert t1 87.6 7.62 12 8.70
## 5 sin fert t2 87.8 7.42 12 8.45
## 6 sin fert t3 87.7 8.14 12 9.28Resumen estadistico visual
#visualizacion
library(ggplot2)
ggplot(datos2)+
aes(tiempo, rto, fill=treatment)+
geom_boxplot()
#Fertilizar y no fertilizar estan teniendo el mismo efecto
Deteccion de outliers
datos2 %>%
group_by(treatment, tiempo) %>%
identify_outliers(rto)## [1] treatment tiempo id rto is.outlier is.extreme
## <0 rows> (or 0-length row.names)#0 filas, no hay datos atipicos
Normalidad
datos2 %>%
group_by(treatment, tiempo) %>%
shapiro_test(rto)## # A tibble: 6 × 5
## treatment tiempo variable statistic p
## <fct> <chr> <chr> <dbl> <dbl>
## 1 con fert t1 rto 0.828 0.0200
## 2 con fert t2 rto 0.868 0.0618
## 3 con fert t3 rto 0.887 0.107
## 4 sin fert t1 rto 0.919 0.279
## 5 sin fert t2 rto 0.923 0.316
## 6 sin fert t3 rto 0.886 0.104#Nota: La normalidad se hace sobre los residuales, en este caso se esta haciendo sobre los datos de rendimiento ##t1 con fertilizante no presenta datos normales
res.aov <- anova_test(
data = datos2,
dv = rto,
wid = id,
within = c(treatment,
tiempo)
)
get_anova_table(res.aov)## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 treatment 1.00 11.00 15.541 2.00e-03 * 0.059
## 2 tiempo 1.31 14.37 27.369 5.03e-05 * 0.049
## 3 treatment:tiempo 2.00 22.00 30.424 4.63e-07 * 0.050#Nota: lo primero que se debe interpretar es la interaccion (tratamiento:tiempo) #p-valor interaccion 4.63e-07. Como es menor al 5% significa que SI hay interaccion. Cuando SI hay interaccion no se ve los otros factores (No se hacen comparaciones)
#grafico de interaccion
interaction.plot(datos2$tiempo,
datos2$treatment,
datos2$rto)
datos2 %>%
group_by(tiempo, treatment) %>%
summarise(mean_rto = mean(rto)) %>%
ggplot()+
aes(tiempo, mean_rto,
color=treatment,
group=treatment)+
geom_point(size=5)+
geom_line(linewidth=3)## `summarise()` has grouped output by 'tiempo'. You can override using the
## `.groups` argument.
*********************************************************************+
Medidas repetidas de 3 vias 09-06-23
#library(rstatix)
#library(datarium)
data("weightloss", package = "datarium")
set.seed(123)
data("weightloss", package = "datarium")
datos3 = weightloss %>%
rename(FQ = diet,
FO = exercises) %>%
pivot_longer(c(t1,t2,t3),
names_to = 'tiempo',
values_to = 'rto')
datos3## # A tibble: 144 × 5
## id FQ FO tiempo rto
## <fct> <fct> <fct> <chr> <dbl>
## 1 1 no no t1 10.4
## 2 1 no no t2 13.2
## 3 1 no no t3 11.6
## 4 2 no no t1 11.6
## 5 2 no no t2 10.7
## 6 2 no no t3 13.2
## 7 3 no no t1 11.4
## 8 3 no no t2 11.1
## 9 3 no no t3 11.4
## 10 4 no no t1 11.1
## # ℹ 134 more rowstable(datos3$FQ, datos3$FO, datos3$tiempo)## , , = t1
##
##
## no yes
## no 12 12
## yes 12 12
##
## , , = t2
##
##
## no yes
## no 12 12
## yes 12 12
##
## , , = t3
##
##
## no yes
## no 12 12
## yes 12 12head(datos3)## # A tibble: 6 × 5
## id FQ FO tiempo rto
## <fct> <fct> <fct> <chr> <dbl>
## 1 1 no no t1 10.4
## 2 1 no no t2 13.2
## 3 1 no no t3 11.6
## 4 2 no no t1 11.6
## 5 2 no no t2 10.7
## 6 2 no no t3 13.2Descriptivo tabular (visualizacion)
datos3 %>% #agrupar por fertilizante quimico, organico y tiempo
group_by(FQ, FO, tiempo) %>% #una vez que los agrupe obtenga un resumen estadistico por grupo
get_summary_stats(rto, type = "mean_sd")## # A tibble: 12 × 7
## FQ FO tiempo variable n mean sd
## <fct> <fct> <chr> <fct> <dbl> <dbl> <dbl>
## 1 no no t1 rto 12 10.9 0.868
## 2 no no t2 rto 12 11.6 1.30
## 3 no no t3 rto 12 11.4 0.935
## 4 no yes t1 rto 12 10.8 1.27
## 5 no yes t2 rto 12 13.4 1.01
## 6 no yes t3 rto 12 16.8 1.53
## 7 yes no t1 rto 12 11.7 0.938
## 8 yes no t2 rto 12 12.4 1.42
## 9 yes no t3 rto 12 13.8 1.43
## 10 yes yes t1 rto 12 11.4 1.09
## 11 yes yes t2 rto 12 13.2 1.22
## 12 yes yes t3 rto 12 14.7 0.625Descriptivo Grafico
library(ggpubr)
ggboxplot(
datos3, x = "FQ", y = "rto",
color = "tiempo", palette = "jco",
facet.by = "FO", short.panel.labs = FALSE
)
#hay atipicos?
datos3 %>%
group_by(FQ, FO, tiempo) %>%
identify_outliers(rto)#se indentifican los datos atipicos del rendimiento ## # A tibble: 5 × 7
## FQ FO tiempo id rto is.outlier is.extreme
## <fct> <fct> <chr> <fct> <dbl> <lgl> <lgl>
## 1 no no t3 2 13.2 TRUE FALSE
## 2 yes no t1 1 10.2 TRUE FALSE
## 3 yes no t1 3 13.2 TRUE FALSE
## 4 yes no t1 4 10.2 TRUE FALSE
## 5 yes no t2 10 15.3 TRUE FALSE#hay datos sospechosos de atipicos, sin embargo is.extreme confirma que no esisten datos atipicos ##NOTA:en el caso: si existiese datos atipicos solo poner una nota que diga “este datos podria afectar los datos de todos los resultados”
Normalidad
ggplot(datos3)+
aes(rto, fill=tiempo)+
geom_density(alpha=0.5)#aqui esta agrupado por tiempo
#Normalidad global
ggplot(datos3)+
aes(rto)+
geom_density()
#pruebas de normalidad (se hace al rendimiento)
shapiro.test(datos3$rto)##
## Shapiro-Wilk normality test
##
## data: datos3$rto
## W = 0.96377, p-value = 0.0007474library(nortest)
lillie.test(datos3$rto)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: datos3$rto
## D = 0.10013, p-value = 0.001229cvm.test(datos3$rto)##
## Cramer-von Mises normality test
##
## data: datos3$rto
## W = 0.23031, p-value = 0.002205sf.test(datos3$rto)##
## Shapiro-Francia normality test
##
## data: datos3$rto
## W = 0.96406, p-value = 0.001295pearson.test(datos3$rto) ##
## Pearson chi-square normality test
##
## data: datos3$rto
## P = 28.5, p-value = 0.004672#Prueba de simetria
library(lawstat)
symmetry.test(datos3$rto)##
## m-out-of-n bootstrap symmetry test by Miao, Gel, and Gastwirth (2006)
##
## data: datos3$rto
## Test statistic = 2.2399, p-value = 0.032
## alternative hypothesis: the distribution is asymmetric.
## sample estimates:
## bootstrap optimal m
## 19#estos datos son simetricos al 1% (3%)
ANALISIS DE VARIANZA
#esta pruab es mas recomendable
res.aov = anova_test(
data = datos3, dv =rto, wid = id,
within = c(FQ, FO, tiempo)
)
get_anova_table(res.aov)## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 FQ 1.00 11.00 6.021 3.20e-02 * 0.028
## 2 FO 1.00 11.00 58.928 9.65e-06 * 0.284
## 3 tiempo 2.00 22.00 110.942 3.22e-12 * 0.541
## 4 FQ:FO 1.00 11.00 75.356 2.98e-06 * 0.157
## 5 FQ:tiempo 1.38 15.17 0.603 5.01e-01 0.013
## 6 FO:tiempo 2.00 22.00 20.826 8.41e-06 * 0.274
## 7 FQ:FO:tiempo 2.00 22.00 14.246 1.07e-04 * 0.147#Hay triple interaccion (1.07e-04 ), ya no se ve ningun otro efecto
interaction.plot(datos3$tiempo,
datos3$FO,
datos3$rto)
interaction.plot(datos3$tiempo,
datos3$FQ,
datos3$rto)
#Interaccion con tabla de medias
datos3 %>%
group_by(tiempo, FQ, FO) %>%
summarise(m_rto = mean(rto)) #tabular completa (es mejor la tabular cruzada(No usar esta tabla))## `summarise()` has grouped output by 'tiempo', 'FQ'. You can override using the
## `.groups` argument.## # A tibble: 12 × 4
## # Groups: tiempo, FQ [6]
## tiempo FQ FO m_rto
## <chr> <fct> <fct> <dbl>
## 1 t1 no no 10.9
## 2 t1 no yes 10.8
## 3 t1 yes no 11.7
## 4 t1 yes yes 11.4
## 5 t2 no no 11.6
## 6 t2 no yes 13.4
## 7 t2 yes no 12.4
## 8 t2 yes yes 13.2
## 9 t3 no no 11.4
## 10 t3 no yes 16.8
## 11 t3 yes no 13.8
## 12 t3 yes yes 14.7#Interaccion con tabla de medias (tabular cruzada)
tbl = tapply(datos3$rto, list(datos3$FQ,
datos3$FO,
datos3$tiempo),
mean)
addmargins(tbl, FUN = mean)## Margins computed over dimensions
## in the following order:
## 1:
## 2:
## 3:## , , t1
##
## no yes mean
## no 10.90917 10.79417 10.85167
## yes 11.74250 11.39333 11.56792
## mean 11.32583 11.09375 11.20979
##
## , , t2
##
## no yes mean
## no 11.56583 13.42083 12.49333
## yes 12.41583 13.22500 12.82042
## mean 11.99083 13.32292 12.65688
##
## , , t3
##
## no yes mean
## no 11.45000 16.81750 14.13375
## yes 13.78667 14.65500 14.22083
## mean 12.61833 15.73625 14.17729
##
## , , mean
##
## no yes mean
## no 11.30833 13.67750 12.49292
## yes 12.64833 13.09111 12.86972
## mean 11.97833 13.38431 12.68132#Organico #El tiempo 2 y 3 estan causando la interaccion debido a que cuaando se une la columna y fila de los margenes (mayor valor de cada margen) el resultado no es el valor mas grande