我想写一个函数批量做ANOVA,但是有一个问题我解决不了,问题出在变量引用上,整个代码写得很正确,因为一切都是“靠脚”计算的,但是当我试图把这段代码插入到函数中时,却不起作用,能不能有人看一下,指出问题出在哪里?
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)
在下面的代码中,一切都运行良好...
# 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")
此代码的结果是
下面的图表
但是这个函数不起作用...对“deparse(substitute(x)”的修复只是部分解决了这个问题,我得到了一个警告:
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.
下面是该函数的完整代码:
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)
我有一个巨大的要求提示如何解决这个问题。
1条答案
按热度按时间lnlaulya1#
问题出在公式上,请尝试添加以下代码: