R/manova-base.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)
library(ggpubr)
library(rstatix)
library(ggplot2)
# One-way repeated measures ANOVA
data("selfesteem", package = "datarium")
head(selfesteem, 10)
#
selfesteem <- selfesteem %>%
gather(key = "time", value = "score", t1, t2, t3) %>%
   convert_as_factor(id, time)
head(selfesteem, 10)
view(selfesteem)
# Summary statistics
selfesteem %>%
   group_by(time) %>%
   get_summary_stats(score, type = "mean_sd")
# Visualization
bxp <- ggboxplot(selfesteem, x = "time", y = "score", add = "point")
bxp
# Check assumptions
# Outliers
selfesteem %>%
   group_by(time) %>%
   identify_outliers(score)
# Normality assumption
selfesteem %>%
   group_by(time) %>%
   shapiro_test(score)
#
ggqqplot(selfesteem, "score", facet.by = "time")
#
# Assumption of sphericity
head(selfesteem)
view(selfesteem)
res.aov <- anova_test(data = selfesteem, dv = score, wid = id, within = time)
get_anova_table(res.aov)
#
# Post-hoc tests: bonferroni
# pairwise comparisons
pwc <- selfesteem %>%
   pairwise_t_test(
      score ~ time, paired = TRUE,
      p.adjust.method = "bonferroni"
   )
pwc
# Report
# Visualization: box plots with p-values
pwc <- pwc %>% add_xy_position(x = "time")
bxp +
   stat_pvalue_manual(pwc) +
   labs(
      subtitle = get_test_label(res.aov, detailed = TRUE),
      caption = get_pwc_label(pwc)
   )
#
# 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)
view(selfesteem2)
#
# 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
# Chek 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)
view(weightloss)
# 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)
   )