Problem with reference to variable in R function-CodePudding

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?

df <- data.frame(
  stringsAsFactors = FALSE,
  id = c(1L,2L,3L,4L, 5L,1L,2L,3L,4L,5L,1L,2L,3L,4L,5L,1L,2L,3L,4L,5L),
  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

But the function doesn't work...

my_function <- function(df, x, y) { 
  # 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"))
  
### The problem arises in this part of the code ###
  # Pairwise comparisons
  pwc <- df %>%
   pairwise_t_test({{y}} ~ {{x}}, paired = TRUE, conf.level = 0.95,
   detailed = TRUE, p.adjust.method = "bonferroni")
  
  output <- list(a,b,pwc)
  return(output)
}

my_function(df, treat, vo2)

CodePudding user response：

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

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

  # Pairwise comparisons
  pwc <- pairwise_t_test(data=df, formula, paired = TRUE, conf.level = 0.95,
                    detailed = TRUE, p.adjust.method = "bonferroni")

CodePudding user response：

Thank you for the tip. This partially solves my problem, although I get a warning::

1: In mean.default(~"vo2", na.rm = TRUE) :  argument nie jest wartością liczbową ani logiczną: zwracanie wartości NA
2: In var(if (is.vector(x) || is.factor(x)) x else as.double(x), na.rm = na.rm) :  NAs introduced by coercion.

Unfortunately, in the further part of my function I wanted to build a graph on this data, but there are also errors related to the reference to variables in the function. I am attaching the whole code, maybe someone can help me solve my problem.

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)

  # 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)   # poziome kreski na koncu wąsów
       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")
    
  output <- list(a,b,pwc2,d)
  return(output)
}

my_function(df, treat, vo2)