如何在 R 中绘制带条件 p 值和置信区间的方差分析图?
How to graph an ANOVA w/ conditional p-values and confidence intervals in R?
我有 df1:
Rate Dogs MHI_2018 Points Level AGE65_MORE P_Elderly
1 0.10791173 0.00000000 59338 236.4064 C 8653 15.56267
2 0.06880040 0.00000000 57588 229.4343 C 44571 20.44335
3 0.08644537 0.00000000 50412 200.8446 C 10548 18.23651
4 0.29591635 0.00000000 29267 116.6016 A 1661 16.38390
5 0.05081301 0.00000000 37365 148.8645 B 3995 20.29980
6 0.02625200 0.00000000 45400 180.8765 D 20247 17.71748
7 0.80321285 0.02974862 39917 159.0319 D 6562 19.52105
8 0.07682852 0.00000000 42132 167.8566 D 5980 22.97173
9 0.18118814 0.00000000 47547 189.4303 B 7411 16.78482
10 0.07787555 0.00000000 39907 158.9920 B 2953 22.99665
11 0.15065913 0.00000000 39201 156.1793 C 2751 20.72316
12 0.33362247 0.00000000 46495 185.2390 B 2915 19.45019
13 0.03652168 0.00000000 49055 195.4382 B 10914 19.92988
14 0.27998133 0.00000000 42423 169.0159 A 2481 23.15446
15 0.05407451 0.00000000 40203 160.1713 A 7790 21.06202
16 0.07233796 0.00000000 39057 155.6056 A 2629 19.01765
17 0.08389061 0.00000000 45796 182.4542 B 15446 18.51106
18 0.05220569 0.00000000 34035 135.5976 B 6921 18.06578
19 0.05603418 0.00000000 39491 157.3347 B 12322 17.26133
20 0.15875536 0.00000000 60367 240.5060 C 12400 15.14282
有
AOV <- aov(Rate~Level, data = df)
TukeyHSD(AOV)
$Level
diff lwr upr p adj
B-A -0.066558621 -0.3783957 0.2452784 0.9272012
C-A -0.061063140 -0.4026635 0.2805372 0.9551663
D-A 0.126520253 -0.2624089 0.5154494 0.7890519
C-B 0.005495482 -0.2848090 0.2958000 0.9999404
D-B 0.193078874 -0.1516699 0.5378277 0.4049948
D-C 0.187583392 -0.1843040 0.5594708 0.4923479
我现在想用均值和置信区间绘制这些数据的图。如果 p_adj < 0.50,我还想在变量之间绘制它。输出看起来像:
一种解决方案是使用 ggsignif
包,但首先您需要准备 TukeyHSD
的输出以供在 ggsignif
中使用:
AOV <- aov(Rate~Level, data = df)
t <-as.data.frame(TukeyHSD(AOV)$Level)
library(tidyverse)
MAX <-df %>% group_by(Level) %>% summarise(Max = max(Rate))
T1 <- t %>% rownames_to_column("Group") %>%
mutate(Start = sub("^(.).*","\1",Group),
End = sub(".*(.)$","\1",Group)) %>%
left_join(.,MAX, by = c("Start" = "Level")) %>%
left_join(.,MAX, by = c("End" = "Level")) %>%
mutate(End = factor(End)) %>%
rowwise() %>%mutate(ypos = max(Max.x, Max.y)*(1+0.25*as.numeric(End)))
Source: local data frame [6 x 10]
Groups: <by row>
# A tibble: 6 x 10
Group diff lwr upr `p adj` Start End Max.x Max.y ypos
<chr> <dbl> <dbl> <dbl> <dbl> <chr> <fct> <dbl> <dbl> <dbl>
1 B-A -0.0666 -0.378 0.245 0.927 B A 0.334 0.296 0.417
2 C-A -0.0611 -0.403 0.281 0.955 C A 0.159 0.296 0.370
3 D-A 0.127 -0.262 0.515 0.789 D A 0.803 0.296 1.00
4 C-B 0.00550 -0.285 0.296 1.00 C B 0.159 0.334 0.500
5 D-B 0.193 -0.152 0.538 0.405 D B 0.803 0.334 1.20
6 D-C 0.188 -0.184 0.559 0.492 D C 0.803 0.159 1.41
现在,您可以绘制数据并根据 T1 数据集添加显着性:
library(ggsignif)
library(ggplot2)
ggplot(df, aes(x = Level, y = Rate))+
geom_jitter(width = 0.2)+
stat_summary(fun.data = "mean_cl_normal", geom = "errorbar", width = 0, color = "red") +
stat_summary(fun = "mean", geom = "errorbar", aes(ymax = ..y.., ymin = ..y..), col = "red", width = 0.5) +
geom_signif(data = subset(T1,`p adj` <0.5), manual = TRUE,
aes(xmax = End, xmin = Start, y_position= ypos, annotations = round(`p adj`,3)))
我有 df1:
Rate Dogs MHI_2018 Points Level AGE65_MORE P_Elderly
1 0.10791173 0.00000000 59338 236.4064 C 8653 15.56267
2 0.06880040 0.00000000 57588 229.4343 C 44571 20.44335
3 0.08644537 0.00000000 50412 200.8446 C 10548 18.23651
4 0.29591635 0.00000000 29267 116.6016 A 1661 16.38390
5 0.05081301 0.00000000 37365 148.8645 B 3995 20.29980
6 0.02625200 0.00000000 45400 180.8765 D 20247 17.71748
7 0.80321285 0.02974862 39917 159.0319 D 6562 19.52105
8 0.07682852 0.00000000 42132 167.8566 D 5980 22.97173
9 0.18118814 0.00000000 47547 189.4303 B 7411 16.78482
10 0.07787555 0.00000000 39907 158.9920 B 2953 22.99665
11 0.15065913 0.00000000 39201 156.1793 C 2751 20.72316
12 0.33362247 0.00000000 46495 185.2390 B 2915 19.45019
13 0.03652168 0.00000000 49055 195.4382 B 10914 19.92988
14 0.27998133 0.00000000 42423 169.0159 A 2481 23.15446
15 0.05407451 0.00000000 40203 160.1713 A 7790 21.06202
16 0.07233796 0.00000000 39057 155.6056 A 2629 19.01765
17 0.08389061 0.00000000 45796 182.4542 B 15446 18.51106
18 0.05220569 0.00000000 34035 135.5976 B 6921 18.06578
19 0.05603418 0.00000000 39491 157.3347 B 12322 17.26133
20 0.15875536 0.00000000 60367 240.5060 C 12400 15.14282
有
AOV <- aov(Rate~Level, data = df)
TukeyHSD(AOV)
$Level
diff lwr upr p adj
B-A -0.066558621 -0.3783957 0.2452784 0.9272012
C-A -0.061063140 -0.4026635 0.2805372 0.9551663
D-A 0.126520253 -0.2624089 0.5154494 0.7890519
C-B 0.005495482 -0.2848090 0.2958000 0.9999404
D-B 0.193078874 -0.1516699 0.5378277 0.4049948
D-C 0.187583392 -0.1843040 0.5594708 0.4923479
我现在想用均值和置信区间绘制这些数据的图。如果 p_adj < 0.50,我还想在变量之间绘制它。输出看起来像:
一种解决方案是使用 ggsignif
包,但首先您需要准备 TukeyHSD
的输出以供在 ggsignif
中使用:
AOV <- aov(Rate~Level, data = df)
t <-as.data.frame(TukeyHSD(AOV)$Level)
library(tidyverse)
MAX <-df %>% group_by(Level) %>% summarise(Max = max(Rate))
T1 <- t %>% rownames_to_column("Group") %>%
mutate(Start = sub("^(.).*","\1",Group),
End = sub(".*(.)$","\1",Group)) %>%
left_join(.,MAX, by = c("Start" = "Level")) %>%
left_join(.,MAX, by = c("End" = "Level")) %>%
mutate(End = factor(End)) %>%
rowwise() %>%mutate(ypos = max(Max.x, Max.y)*(1+0.25*as.numeric(End)))
Source: local data frame [6 x 10]
Groups: <by row>
# A tibble: 6 x 10
Group diff lwr upr `p adj` Start End Max.x Max.y ypos
<chr> <dbl> <dbl> <dbl> <dbl> <chr> <fct> <dbl> <dbl> <dbl>
1 B-A -0.0666 -0.378 0.245 0.927 B A 0.334 0.296 0.417
2 C-A -0.0611 -0.403 0.281 0.955 C A 0.159 0.296 0.370
3 D-A 0.127 -0.262 0.515 0.789 D A 0.803 0.296 1.00
4 C-B 0.00550 -0.285 0.296 1.00 C B 0.159 0.334 0.500
5 D-B 0.193 -0.152 0.538 0.405 D B 0.803 0.334 1.20
6 D-C 0.188 -0.184 0.559 0.492 D C 0.803 0.159 1.41
现在,您可以绘制数据并根据 T1 数据集添加显着性:
library(ggsignif)
library(ggplot2)
ggplot(df, aes(x = Level, y = Rate))+
geom_jitter(width = 0.2)+
stat_summary(fun.data = "mean_cl_normal", geom = "errorbar", width = 0, color = "red") +
stat_summary(fun = "mean", geom = "errorbar", aes(ymax = ..y.., ymin = ..y..), col = "red", width = 0.5) +
geom_signif(data = subset(T1,`p adj` <0.5), manual = TRUE,
aes(xmax = End, xmin = Start, y_position= ypos, annotations = round(`p adj`,3)))