R/manova-RedRio.R

In carlosperezoft/hipervizr: Dashboard para Hipercartas y Cartas de Control de Red Rio.

# autor --------------------------------------
# carlos.perez7@udea.edu.co
# carlos.perezoft@gmail.com
# 31/03/2022 4:40:35 p. m.
# autor --------------------------------------
#
library(tidyverse) # v1.3.1
library(ggpubr)    # v0.4.0
library(rstatix)   # v0.7.0
library(ggplot2)   # v3.3.5
#
library(readxl, quietly=TRUE) # v1.3.1
# Ruta base para los archivos de datos de la app HIPERVIZ:
hiperviz_data_path <- "C:\\Temp\\"
#
#
mrE1Data <- read_excel(paste0(hiperviz_data_path,"CONSOLIDADO_E1_CONV_825_c15mins.xlsx"))
head(mrE1Data, 10)
dbRepetidas <- mrE1Data[c("t", "DIA_MES","MES","ANIO","PH")]
#
dbRepetidas <- dbRepetidas %>% convert_as_factor(t,MES,DIA_MES)
#
# View funcion que carga una ventana tipo "data table" con los datos indicados:
head(dbRepetidas, 10)
view(dbRepetidas)
#
# One-way repeated measures ANOVA
statRep <- dbRepetidas %>%
            group_by(MES,DIA_MES,ANIO) %>%
            get_summary_stats(PH, type = "five_number")
##
view(statRep)
#
boxRep <- ggboxplot(dbRepetidas, x = "MES", y = "PH", add = "point", color = "MES", palette = "d3")
boxRep
#
# Outliers: Se puede agrupar por la combinacion MES, DIA_MES o individual, segun el analisis:
atipicosRep <- dbRepetidas %>%
               group_by(MES, ANIO) %>%
               identify_outliers(PH)
view(atipicosRep)
#

###
# Normality assumption
# Outliers: Se puede agrupar por la combinacion MES, DIA_MES o individual, segun el analisis:
normalRep <- dbRepetidas %>%
         group_by(MES,ANIO) %>%
         shapiro_test(PH)
view(normalRep)
#
ggqqplot(dbRepetidas, "PH", facet.by = "MES")
ggqqplot(dbRepetidas, x = "PH", conf.int = TRUE, facet.by = "MES", color = "MES", palette = "d3", shape = "o") # palette = "tron/simpsons/futurama/startrek"
#
# Assumption of sphericity
head(dbRepetidas, 12)
# Se especifica como: within la combinacion de MES+DIA_MES ya que da una combinacion unica con el id "t".
# Se debe usar c(MES,DIA_MES), para que MES represente el "tratamiento":
redrio.aov <- anova_test(data = dbRepetidas, dv = PH, wid = t, within = c(MES,DIA_MES))
anovaTable <- get_anova_table(redrio.aov) # Calcula la prueba de esfericidad!!
view(anovaTable)
#
# Post-hoc tests: bonferroni
# pairwise comparisons
pwt.rep <- dbRepetidas %>% group_by(DIA_MES) %>%
   pairwise_t_test(
      PH ~ MES, paired = TRUE, #DIA_MES
      p.adjust.method = "bonferroni"
   )
view(pwt.rep)
# Report
# Visualization: box plots with p-values
pwt.rep <- pwt.rep %>% add_xy_position(x="DIA_MES") #
boxRep +
   stat_pvalue_manual(pwt.rep) +
   labs(
      caption = get_pwc_label(pwt.rep)
   )
#
# Two-way repeated measures ANOVA
#
# The two-way repeated measures ANOVA can be performed in order to determine
# whether there is a significant interaction between diet and time on the self-esteem score.
#
# Wide format: Data preparation
set.seed(123)
data("selfesteem2", package = "datarium")
selfesteem2 %>% sample_n_by(treatment, size = 1)
#
# Gather the columns t1, t2 and t3 into long format.
# Convert id and time into factor variables
selfesteem2 <- selfesteem2 %>%
   gather(key = "time", value = "score", t1, t2, t3) %>%
   convert_as_factor(id, time)
# Inspect some random rows of the data by groups
set.seed(123)
selfesteem2 %>% sample_n_by(treatment, time, size = 1)
#
# Summary statistics
selfesteem2 %>%
   group_by(treatment, time) %>%
   get_summary_stats(score, type = "mean_sd")
# Visualization
bxp <- ggboxplot(
   selfesteem2, x = "time", y = "score",
   color = "treatment", palette = "jco"
)
bxp
# Check assumptions: Outliers
selfesteem2 %>%
   group_by(treatment, time) %>%
   identify_outliers(score)
# Normality assumption
selfesteem2 %>%
   group_by(treatment, time) %>%
   shapiro_test(score)
# Create QQ plot for each cell of design:
ggqqplot(selfesteem2, "score", ggtheme = theme_bw()) +
   facet_grid(time ~ treatment, labeller = "label_both")
# Computation: ANOVA test:
res.aov <- anova_test(
   data = selfesteem2, dv = score, wid = id,
   within = c(treatment, time)
)
get_anova_table(res.aov)
#
# Post-hoc tests
# Procedure for a significant two-way interaction
# Effect of treatment at each time point
one.way <- selfesteem2 %>%
   group_by(time) %>%
   anova_test(dv = score, wid = id, within = treatment) %>%
   get_anova_table() %>%
   adjust_pvalue(method = "bonferroni")
one.way
# Pairwise comparisons between treatment groups
pwc <- selfesteem2 %>%
   group_by(time) %>%
   pairwise_t_test(
      score ~ treatment, paired = TRUE,
      p.adjust.method = "bonferroni"
   )
pwc
# Effect of time at each level of treatment
one.way2 <- selfesteem2 %>%
   group_by(treatment) %>%
   anova_test(dv = score, wid = id, within = time) %>%
   get_anova_table() %>%
   adjust_pvalue(method = "bonferroni")
# Pairwise comparisons between time points
pwc2 <- selfesteem2 %>%
   group_by(treatment) %>%
   pairwise_t_test(
      score ~ time, paired = TRUE,
      p.adjust.method = "bonferroni"
   )
pwc2
# Procedure for non-significant two-way interaction
# comparisons for treatment variable
selfesteem2 %>%
   pairwise_t_test(
      score ~ treatment, paired = TRUE,
      p.adjust.method = "bonferroni"
   )
# comparisons for time variable
selfesteem2 %>%
   pairwise_t_test(
      score ~ time, paired = TRUE,
      p.adjust.method = "bonferroni"
   )
# Report:
# Visualization: box plots with p-values
pwc <- pwc %>% add_xy_position(x = "time")
bxp +
   stat_pvalue_manual(pwc, tip.length = 0, hide.ns = TRUE) +
   labs(
      subtitle = get_test_label(res.aov, detailed = TRUE),
      caption = get_pwc_label(pwc)
   )
#
# Three-way repeated measures ANOVA
#
# Data preparation:
# Wide format
set.seed(123)
data("weightloss", package = "datarium")
weightloss %>% sample_n_by(diet, exercises, size = 1)
#
# Gather the columns t1, t2 and t3 into long format.
# Convert id and time into factor variables
weightloss <- weightloss %>%
   gather(key = "time", value = "score", t1, t2, t3) %>%
   convert_as_factor(id, time)
# Inspect some random rows of the data by groups
set.seed(123)
weightloss %>% sample_n_by(diet, exercises, time, size = 1)
# Summary statistics
weightloss %>%
   group_by(diet, exercises, time) %>%
   get_summary_stats(score, type = "mean_sd")
# Visualization:
bxp <- ggboxplot(
   weightloss, x = "exercises", y = "score",
   color = "time", palette = "jco",
   facet.by = "diet", short.panel.labs = FALSE
)
bxp
# Check assumptions: outliers
weightloss %>%
   group_by(diet, exercises, time) %>%
   identify_outliers(score)
# Normality assumption
weightloss %>%
   group_by(diet, exercises, time) %>%
   shapiro_test(score)
# QQ plot:
ggqqplot(weightloss, "score", ggtheme = theme_bw()) +
   facet_grid(diet + exercises ~ time, labeller = "label_both")
#Computation
res.aov <- anova_test(
   data = weightloss, dv = score, wid = id,
   within = c(diet, exercises, time)
)
get_anova_table(res.aov)
# Post-hoc tests:
# Compute simple two-way interaction
# Two-way ANOVA at each diet level
two.way <- weightloss %>%
   group_by(diet) %>%
   anova_test(dv = score, wid = id, within = c(exercises, time))
two.way
# Extract anova table
get_anova_table(two.way)
#Compute simple simple main effect
# Effect of time at each diet X exercises cells
time.effect <- weightloss %>%
   group_by(diet, exercises) %>%
   anova_test(dv = score, wid = id, within = time)
time.effect
# Extract anova table
get_anova_table(time.effect) %>%
   filter(diet == "no")
# Compute simple simple comparisons
# Pairwise comparisons
pwc <- weightloss %>%
   group_by(diet, exercises) %>%
   pairwise_t_test(score ~ time, paired = TRUE, p.adjust.method = "bonferroni") %>%
   select(-df, -statistic) # Remove details
# Show comparison results for "diet:no,exercises:yes" groups
pwc %>% filter(diet == "no", exercises == "yes") %>%
   select(-p)     # remove p columns
# Report
# Visualization: box plots with p-values
pwc <- pwc %>% add_xy_position(x = "exercises")
pwc.filtered <- pwc %>%
   filter(diet == "no", exercises == "yes")
bxp +
   stat_pvalue_manual(pwc.filtered, tip.length = 0, hide.ns = TRUE) +
   labs(
      subtitle = get_test_label(res.aov, detailed = TRUE),
      caption = get_pwc_label(pwc)
   )