0

I'd like to write a function to do ANOVA in batch, but there's a problem I can't solve. The problem is with the variable reference. The whole code is written correctly, because everything is calculated "on foot", but when I try to insert this code into the function, it doesn't work. Can someone take a look and point out where the problem is?

library(tidyverse)
library(ggpubr)
library(rstatix)
library(ggprism)

df <- data.frame(
  stringsAsFactors = FALSE,
  id = c(1,2,3,4,5,1,2,3,4,5,1,2,3,4,5,1,2,3,4,5),
  treat = c("o","o","o","o","o","j","j","j","j","j","z","z","z","z","z","w","w","w","w","w"),
  vo2 = c("47.48","42.74","45.23","51.65","49.11","51.00","43.82","49.88","54.61","52.20","51.31",
          "47.56","50.69","54.88","55.01","51.89","46.10","50.98","53.62","52.77"))

df$vo2 <- as.numeric(df$vo2)
df$treat <- factor(df$treat)

Everything works fine in the code below...

# Summary
group_by(df, treat) %>% 
  summarise(
    N = n(),
    Mean = mean(vo2, na.rm = TRUE),
    Sd = sd(vo2, na.rm = TRUE))

# ANOVA
res.aov <- anova_test(dv = vo2, wid = id, within = treat, data = df)
get_anova_table(res.aov, correction = c("auto"))

# Pairwise comparisons
pwc <- df %>%
  pairwise_t_test(vo2 ~ treat, paired = TRUE, conf.level = 0.95,
                  detailed = TRUE, p.adjust.method = "bonferroni")
pwc

pwc <- pwc %>% add_xy_position(x = "treat")

ggplot(df, aes(x = treat, y = vo2)) +
  stat_boxplot(aes(x = treat, y = vo2, color = treat), geom = 'errorbar', coef=1.5, width=0.4, linetype = 1) +
  geom_boxplot(aes(x = treat, y = vo2, color = treat, fill = treat)) +
  geom_jitter(aes(x = treat, y = vo2, color = treat, fill = treat), width = 0.2) +
  stat_summary(fun = mean, geom = "point", shape = 0, size = 2, color = "black", stroke = 1) +
  #xlab(deparse(substitute(x))) +   ylab(deparse(substitute(y))) +
  theme_prism(base_size = 14) + scale_x_discrete(guide = "prism_bracket") +
  scale_fill_prism(palette = "floral") +   scale_colour_prism(palette = "floral") +
  scale_y_continuous(expand = expansion(mult = c(0.02, 0.05))) +
  stat_pvalue_manual(pwc, tip.length = 0, hide.ns = TRUE) +
  labs(
    subtitle = get_test_label(res.aov, detailed = TRUE),
    caption = get_pwc_label(pwc)) +
  theme(legend.position = "NULL")


The result of this code is the chart below enter image description here

But the function doesn't work... The fix for "deparse(substitute(x)" only partially solved the problem, and I get a warning:

1: In mean.default(~"vo2", na.rm = TRUE): argument is not numeric or logical: returning an
2: In var(if (is.vector(x) || is.factor(x)) x else as.double(x), na.rm = na.rm) : NAs introduced by coercion.

Below is the full code for this function:

my_function <- function(df, x, y) { 
  
  x <- deparse(substitute(x))
  y <- deparse(substitute(y))
  
  formula <- as.formula(paste0(y, "~", x))
  
  # Summary
a <- group_by(df, {{x}}) %>% 
    summarise(
      N = n(),
      Mean = mean({{y}}, na.rm = TRUE),
      Sd = sd({{y}}, na.rm = TRUE)
    )
  # ANOVA
  res.aov <<- anova_test(dv = {{y}}, wid = id, within = {{x}}, data = df)
b <- get_anova_table(res.aov, correction = c("auto"))
  
  # Pairwise comparisons
pwc <-  pairwise_t_test(data = df, formula, paired = TRUE, conf.level = 0.95,
                          detailed = TRUE, p.adjust.method = "bonferroni")
  
pwc2 <<- pwc %>% add_xy_position(x = {{x}})
  
  # Plot
d <- ggplot(df, aes(x = {{x}}, y = {{y}})) +
       stat_boxplot(aes(x = {{x}}, y = {{y}}, color = {{x}}), geom = 'errorbar', coef=1.5, width=0.4, linetype = 1) +
       geom_boxplot(aes(x = {{x}}, y = {{y}}, color = {{x}}, fill = {{x}})) +
       geom_jitter(aes(x = {{x}}, y = {{y}}, color = {{x}}, fill = {{x}}), width = 0.2) +
    stat_summary(fun = mean, geom = "point", shape = 0, size = 2, color = "black", stroke = 1) +
    stat_pvalue_manual(pwc2, tip.length = 0, hide.ns = TRUE) +
      xlab(deparse(substitute(x))) + ylab(deparse(substitute(y))) +
      scale_x_discrete(guide = "prism_bracket") +
    theme_prism(base_size = 14) +
    theme(legend.position = "NULL")
  
  #ggsave(paste0(deparse(substitute(x)), "_",
  # deparse(substitute(y)), ".png"), width=160, height=90, units="mm", dpi=600)
  
   output <- list(a,b,pwc2,d)
  return(output)
}

my_function(df, treat, vo2)

I have a huge request for tips on how to solve this.

2
  • See my updated answer. The code should be put immediately before the pairwise_t_test command. Commented Dec 30, 2022 at 2:59
  • @Edward Thank you very much for your help in solving this problem. Now the whole function works perfectly as I wanted it to work :-). In the code for calculating pwc, a small correction should be made, namely replacing {{x}} with {{x1}} - like this: pwc2 <<- pwc %>% add_xy_position(x = {{x1}}). Many thanks again!!! Commented Dec 30, 2022 at 10:18

1 Answer 1

1

The problem lies in the formula. Try adding the following code:

  x1 <- deparse(substitute(x))
  y1 <- deparse(substitute(y))
  
  formula <- as.formula(paste0(y1, "~", x1))

  # Pairwise comparisons
  pwc <- pairwise_t_test(data=df, formula, paired = TRUE, conf.level = 0.95,
                    detailed = TRUE, p.adjust.method = "bonferroni")
Sign up to request clarification or add additional context in comments.

Comments

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.