Anova Medias Repetidas 2 (viabilidad bacteriana)
2023-08-18
1.-INTRODUCCIÓN
En el Anova de medidas repetidas se analizan los resultados obtenidos en un diseño experimental en donde se ha manipulado una única variable independiente (un único factor) con 2 o más niveles pero de forma intra-sujeto. Esto viene a significar, que todos los individuos (o unidades de observación), han pasado por todos los niveles del factor. A este tipo de diseños también se les conoce como diseños de medidas repetidas en el sentido de que a cada sujeto se le repite la medición de la variable dependiente de respuesta en diversas condiciones, tantas como niveles tenga el factor manipulado. También se les conoce como diseños de medidas dependientes debido a que las puntuaciones de un mismo sujeto muestran dependencia estadística entre ellas, están relacionadas.
Se trata de comprobar, por ejemplo, desde un punto de vista de contrastes paramétricos, si las medias de una medición de una variable dependiente continua, son las mismas en distintos momentos del tiempo (3 o más), por ejemplo si un tratamiento de combinación de dieta y ejercicio (las 2 cosas como un todo) tiene la correspondiente contraprestación en pérdida de peso en 4 instantes del tiempo, esto es, se trata de 4 muestras relacionadas, dependientes o pareadas, tomadas sobre los mismos 10 individuos. Es algo parecido a la prueba de la T de Student para muestras relacionadas en un antes y un después, pero esta vez con 3 o más instantes de tiempo, en lugar de 2.
2.- EJERCICIO
2.1.- ANÁLISIS DE LOS DATOS
library(tidyverse)## Warning: package 'tidyverse' was built under R version 4.2.3## Warning: package 'ggplot2' was built under R version 4.2.3## Warning: package 'tibble' was built under R version 4.2.3## Warning: package 'tidyr' was built under R version 4.2.3## Warning: package 'readr' was built under R version 4.2.3## Warning: package 'purrr' was built under R version 4.2.3## Warning: package 'dplyr' was built under R version 4.2.3## Warning: package 'forcats' was built under R version 4.2.3## Warning: package 'lubridate' was built under R version 4.2.3## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.2 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.2 ✔ tibble 3.2.1
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errorslibrary(ggplot2)
library(car)## Warning: package 'car' was built under R version 4.2.3## Loading required package: carData## Warning: package 'carData' was built under R version 4.2.3##
## Attaching package: 'car'
##
## The following object is masked from 'package:dplyr':
##
## recode
##
## The following object is masked from 'package:purrr':
##
## somedf <- read_csv("https://docs.google.com/spreadsheets/d/e/2PACX-1vTGUepm0ddAefpaILx5QR9eCBgvnY-st6FBfVf39W1meHwGGsXvb63iVhNFZ3DvwXCM9ZsFaVo4sdyp/pub?gid=1906853507&single=true&output=csv")## Rows: 126 Columns: 4
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (3): Tipo, Tiempo, Viabilidad bacteriana
## dbl (1): ID
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.cols(
Tipo = col_character(),
Tiempo = col_character(),
`Viabilidad bacteriana` = col_number()
)## cols(
## Tipo = col_character(),
## Tiempo = col_character(),
## `Viabilidad bacteriana` = col_number()
## )str(df)## spc_tbl_ [126 × 4] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ ID : num [1:126] 1 2 3 4 5 6 7 8 9 10 ...
## $ Tipo : chr [1:126] "Placebo" "Placebo" "Placebo" "Placebo" ...
## $ Tiempo : chr [1:126] "T0" "T1" "T2" "T3" ...
## $ Viabilidad bacteriana: chr [1:126] "100" "77.29" "51.36" "51.3" ...
## - attr(*, "spec")=
## .. cols(
## .. ID = col_double(),
## .. Tipo = col_character(),
## .. Tiempo = col_character(),
## .. `Viabilidad bacteriana` = col_character()
## .. )
## - attr(*, "problems")=<externalptr>df## # A tibble: 126 × 4
## ID Tipo Tiempo `Viabilidad bacteriana`
## <dbl> <chr> <chr> <chr>
## 1 1 Placebo T0 100
## 2 2 Placebo T1 77.29
## 3 3 Placebo T2 51.36
## 4 4 Placebo T3 51.3
## 5 5 Placebo T4 89.42
## 6 6 Placebo T5 80.71
## 7 7 Placebo T6 99.69
## 8 8 Placebo T0 100
## 9 9 Placebo T1 99.78
## 10 10 Placebo T2 93.58
## # ℹ 116 more rows2.- OUTLIERS
df$ 'Viabilidad bacteriana' <- as.numeric(df$'Viabilidad bacteriana')str(df)## spc_tbl_ [126 × 4] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ ID : num [1:126] 1 2 3 4 5 6 7 8 9 10 ...
## $ Tipo : chr [1:126] "Placebo" "Placebo" "Placebo" "Placebo" ...
## $ Tiempo : chr [1:126] "T0" "T1" "T2" "T3" ...
## $ Viabilidad bacteriana: num [1:126] 100 77.3 51.4 51.3 89.4 ...
## - attr(*, "spec")=
## .. cols(
## .. ID = col_double(),
## .. Tipo = col_character(),
## .. Tiempo = col_character(),
## .. `Viabilidad bacteriana` = col_character()
## .. )
## - attr(*, "problems")=<externalptr>df %>%
ggplot(aes(x= df$'Viabilidad bacteriana'))+
geom_histogram()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
df %>%
ggplot(aes(x= df$'Viabilidad bacteriana'))+
geom_boxplot()
Se identifican valores extremos.
boxplot.stats(df$'Viabilidad bacteriana')$out## [1] 213.32 -52.77 -72.68 -63.31boxplot.stats(df$'Viabilidad bacteriana')## $stats
## [1] -17.16 44.89 77.96 100.00 161.55
##
## $n
## [1] 126
##
## $conf
## [1] 70.20285 85.71715
##
## $out
## [1] 213.32 -52.77 -72.68 -63.31df %>%
ggplot(aes(x=`Viabilidad bacteriana`))+
geom_histogram()+
facet_grid(Tipo ~ Tiempo)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
summary(df$`Viabilidad bacteriana`)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -72.68 44.96 77.96 70.47 100.00 213.323.- NAS
any(is.na(df$`Viabilidad bacteriana`))## [1] FALSENo tenemos datos ausentes.
4.- DESCRIPTIVOS
df %>%
group_by(Tipo, Tiempo) %>%
summarise(n=n(),
mean = mean(`Viabilidad bacteriana`),
var = var(`Viabilidad bacteriana`),
sd = sd(`Viabilidad bacteriana`
))## `summarise()` has grouped output by 'Tipo'. You can override using the
## `.groups` argument.## # A tibble: 14 × 6
## # Groups: Tipo [2]
## Tipo Tiempo n mean var sd
## <chr> <chr> <int> <dbl> <dbl> <dbl>
## 1 Enjuague T0 13 100 0 0
## 2 Enjuague T1 13 47.2 806. 28.4
## 3 Enjuague T2 13 48.3 2731. 52.3
## 4 Enjuague T3 13 47.9 3014. 54.9
## 5 Enjuague T4 13 61.3 1762. 42.0
## 6 Enjuague T5 13 74.2 2244. 47.4
## 7 Enjuague T6 13 68.9 5109. 71.5
## 8 Placebo T0 5 100 0 0
## 9 Placebo T1 5 82.4 432. 20.8
## 10 Placebo T2 5 80.2 486. 22.0
## 11 Placebo T3 5 79.0 357. 18.9
## 12 Placebo T4 5 74.0 570. 23.9
## 13 Placebo T5 5 97.8 1422. 37.7
## 14 Placebo T6 5 98.2 807. 28.4De forma general, apreciamos diferentes medias y dispersiones, analizando el tratamiento y el efecto temporal.
5.- ANOVA
modelo = aov(df$`Viabilidad bacteriana`~ df$Tipo*df$Tiempo)
summary(modelo)## Df Sum Sq Mean Sq F value Pr(>F)
## df$Tipo 1 13839 13839 7.587 0.00686 **
## df$Tiempo 6 28910 4818 2.642 0.01961 *
## df$Tipo:df$Tiempo 6 3491 582 0.319 0.92591
## Residuals 112 204293 1824
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Tipo: El p valor es 0.00686 < 0.05. Tenemos medias diferentes para cada tipo de tratamiento.
Tiempo: El p valor es 0.019 < 0.05. Tenemos medias diferentes para cada eje temporal.
Interacción: el p valor es 0.32 > 0.05. Las medias son iguales en la interacción entre tiempo y tratamiento.
CONTESTACIÓN A LA PREGUNTA El efecto del tratamiento tiene una interacción no significativa o irrelevante, es decir, que la diferencia bacteriana entre el tratamiento y el placebo ofrece los mismos resultados a lo largo del tiempo. Por ello, el tratamiento no muestra efectos significativos con el paso del tiempo.
plot(modelo)
residuo <- modelo$residualsshapiro.test(residuo)##
## Shapiro-Wilk normality test
##
## data: residuo
## W = 0.96651, p-value = 0.003251El p valor es 0.0032 < 0.05. Tenemos un problema de falta de normalidad en los datos.
df_residuo <- data.frame(residuo=residuo)
# Histograma con ggplot2
p <- ggplot(df_residuo, aes(x=residuo)) +
geom_histogram(aes(y=..density..), bins=30, fill="blue", alpha=0.5) +
geom_density(color="red", size=1.2)## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.# Superponer curva de normalidad
mu <- mean(df_residuo$residuo, na.rm=TRUE)
sigma <- sd(df_residuo$residuo, na.rm=TRUE)
p + stat_function(fun=dnorm, args=list(mean=mu, sd=sigma), color="green", size=1.2) +
labs(title="Histograma con Curva de Normalidad",
x="Residuo", y="Densidad") +
theme_minimal()## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
Se observa que en nuestro histograma, el kernel (la línea roja) muestra
una distribución más apuntada que una distribución normal (kernel
verde). Por tener una distribución con asimetría a la derecha y más
apuntada (leptocúrtica), no se acepta la normalidad.
library(car)df_residuo$grupo <- ifelse(df_residuo$residuo > 0, "Positivo", "Negativo")leveneTest(residuo ~ grupo, data=df_residuo)## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.0209 0.8853
## 124El p valor es 0.88 > 0.05. Tenemos evidencia acerca de la homocedasticidad.
Tenemos un problema de falta de normalidad, pero los resultados del ANOVA son robustos en este problema. Los resultados son valorables, aunque convendría realizar una prueba de Kruskal-Wallis en caso de que los datos fueran reducidos (n < 10 en cada factor).
6.- INFORME DE LA SESIÓN
sesion_info <- devtools::session_info()
dplyr::select(
tibble::as_tibble(sesion_info$packages),
c(package, loadedversion, source)
)## # A tibble: 81 × 3
## package loadedversion source
## <chr> <chr> <chr>
## 1 abind 1.4-5 CRAN (R 4.2.0)
## 2 bit 4.0.5 CRAN (R 4.2.3)
## 3 bit64 4.0.5 CRAN (R 4.2.3)
## 4 bslib 0.4.2 CRAN (R 4.2.2)
## 5 cachem 1.0.6 CRAN (R 4.2.2)
## 6 callr 3.7.3 CRAN (R 4.2.3)
## 7 car 3.1-2 CRAN (R 4.2.3)
## 8 carData 3.0-5 CRAN (R 4.2.3)
## 9 cli 3.6.0 CRAN (R 4.2.2)
## 10 colorspace 2.1-0 CRAN (R 4.2.3)
## # ℹ 71 more rows