Reputation: 89
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
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.
Upvotes: 0
Views: 165
Reputation: 18493
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")
Upvotes: 1