Scripts Overview
Welcome to the phdscripts wiki!
Add p-values to faceted ggplot in R Studio
Problem: stat_compare_means conflicts with ggplot2()'s fill() function in R studio when using facet_wrap. In addition, stat_compare_means statistical tests do not allow for apprehension of p-values obtained from datasets analyzed with higher order models like the glm() Poisson distribution. In addition, the p-values generated from simple t-tests tend to conflict with p-values from an lm().
Purpose: To generate accurate plots with p-values affixed indicating pairwise comparisons.
Method: Adding p-Values Manually to plots in R.
Dependencies: Tidyverse, ggpubr, ggplot2
Step One: Create the data frame listing the column names and values. Because ggpubr() statistical tests are no capable of pulling p-values from higher order statistical tests like glm(Y~X, family=Poisson, data=data), the results from each statistical test will need to be provided into data frame manually.
stat.testAll<- tibble::tribble(~group1, ~group2,~Snp_Type, ~glmModelp, "Middle","Aged","Transversion","3.40e-05", "Middle", "NR","Transversion","3.22e-05", "Middle", "R", "Transversion","<2e-16","R", "Aged","Transversion","1.54e-07", "R","NR","Transversion","5.09e-13","NR","Aged","Transversion","0.577","Middle","Aged","Transition","2.83e-14", "Middle", "NR", "Transition","1.32e-15", "Middle", "R", "Transition", "1.54e-15","R","Aged","Transition","0.891", "R","NR","Transition","0.675","NR","Aged","Transition","0.607")
View dataframe
> stat.testAll
# A tibble: 12 x 4
group1 group2 Snp_Type glmModelp
<chr> <chr> <chr> <chr>
1 Middle Aged Transversion 3.40e-05
2 Middle NR Transversion 3.22e-05
3 Middle R Transversion <2e-16
4 R Aged Transversion 1.54e-07
5 R NR Transversion 5.09e-13
6 NR Aged Transversion 0.577
7 Middle Aged Transition 2.83e-14
8 Middle NR Transition 1.32e-15
9 Middle R Transition 1.54e-15
10 R Aged Transition 0.891
11 R NR Transition 0.675
12 NR Aged Transition 0.607
Step Two: Add a column from which ggplot can use to plot the p-value label onto the final plot using mutate(). The number provided indicate the position along the y-axis the labels will be placed. This provides the "stacked" effect. If there are 2 facets, ensure the values are duplicated.
p1<-stat.testAll %>% mutate(y.position = rep(c(80, 87, 94, 101, 108, 115), 2))
View merged dataframe
> p1
# A tibble: 12 x 5
group1 group2 Snp_Type glmModelp y.position
<chr> <chr> <chr> <chr> <dbl>
1 Middle Aged Transversion 3.40e-05 80
2 Middle NR Transversion 3.22e-05 87
3 Middle R Transversion <2e-16 94
4 R Aged Transversion 1.54e-07 101
5 R NR Transversion 5.09e-13 108
6 NR Aged Transversion 0.577 115
7 Middle Aged Transition 2.83e-14 80
8 Middle NR Transition 1.32e-15 87
9 Middle R Transition 1.54e-15 94
10 R Aged Transition 0.891 101
11 R NR Transition 0.675 108
12 NR Aged Transition 0.607 115
Step Three: Use stat_pvalue_manual() to select the merged dataframe from which ggplot can pull the p-values.
ggplot(snps_brain, aes(x=factor(c(Group),levels=c("NR","R","Middle","Aged")), y=SnpCount, fill=Snp_Type))+ geom_bar(stat="summary", position="dodge",fun="mean")+ labs(title="MtDNA Transitions and Transversions in the Mus Brain\n")+xlab("\nTreatment Group")+ ylab("Average SNP Count")+theme(plot.title = element_text(hjust = 0.5),panel.grid = element_blank(),panel.background = element_blank())+theme(plot.title = element_text(face="bold",size=14, hjust = 0.5),axis.title.y=element_text(size=10, angle=90, face="bold"))+theme_classic()+geom_errorbar(aes(x=Group, ymin=Mean, ymax=Mean+STD_group), width=.2, position=position_dodge(.9))+theme(plot.title = element_text(hjust = 0.5, face="bold"))+theme(axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)))+scale_y_continuous(limits = c(0, 120),expand=c(0,0))+scale_fill_manual("SNP Type",values=c("Transition"="Khaki3","Transversion"="steelblue4"))+facet_wrap(~Snp_Type)+stat_pvalue_manual(p1, label = "glmModelp", tip.length = 0.01)
Resulting plot
