I'm attempting to visualize main effects and interactions from Tukey HSD pairwise comparisons after a 3x12 ANOVA in R. I've been able to add the pairwise comparisons of the second (12-level) factor using code similar to that below. (Below is a much smaller subset of the data to make y'all's life a little easier.)
#Data
df <- data.frame(
ftwo = c("b","a","b","a", "b","a", "b", "a","b","b","a", "a","b","a", "b","a","b","a","b", "b","a","b","a","a","b","b","a","b", "a","a","b","a","b","a", "b","a","a","b","a","a","a","b","b", "a","b", "a","b","a", "b", "a","b","a","b", "a", "b","a","b","a","b","a", "b","b", "a"),
fone=c(1,1,1,1, 1,2, 1,1, 2,1,1,1, 1, 2,1,1, 2,1,2,2,2, 1,1,1,2,2,2,2, 2,1,2,1,1,1, 1,1,2,2, 2,2,1,2, 2,1,2,1, 1,1,2,2,2,2,1,2,1, 2,2,2,2,2,2,1,2),
score=c(1.44337567,0.19409079,0.00000000,0.00000000,0.58227236, -0.97045393,1.19593984,-0.77695466,0.25898489,3.17542648,-0.55277080,-1.57348553, 1.20738671, -0.31773335,-0.65756364,0.00000000,0.00000000,-2.32163735,0.33166248,-0.42817442,-0.42817442,-0.28867514, -0.28867514, -0.13151273,-1.36618661, 0.51176632,0.51176632,-0.28867514,-1.62706279, 0.06043226,1.51080658,-0.28867513, 1.18450885,-1.65831240, 1.65831239,-0.89293744,-1.17260394,0.09255797,-2.12883336,0.00000000,-1.59709882,0.19963735,0.00000000,0.00000000,0.00000000, -0.66687291, 0.99818676,-0.19963735,-0.62279381, 1.08720630,-1.02734847,-0.30216132,-0.85634884, -0.67700320,1.44337567,-0.86602540, 0.71115900, -0.88894875,0.42817442,0.42817442,0.62323196,0.26596035,-1.10183574))
df <- df %>% mutate(
fone = as.factor(fone),
ftwo = as.factor(ftwo)
)
#Pairwise comparisons
df %>% tukey_hsd(score~fone*ftwo)
#Add to plot
tukey <- df %>% tukey_hsd(score~fone*ftwo) %>% filter(term=="ftwo")
df %>% ggplot() +
geom_boxplot(aes(x=ftwo,y=score,fill=fone)) +
stat_pvalue_manual(tukey,label="p.adj",y.position = c(4),size=2)
Example with main effects of the second factor
However, when trying to add the interactions in a similar manner, they don't map onto the plot correctly. If there's some way I should alter the test and/or plot to make them "speak" to each other, or if there's a better approach to this problem, I'd be grateful to hear it.
tukeyfull <- df %>% tukey_hsd(score~fone*ftwo)
df %>% ggplot() +
geom_boxplot(aes(x=ftwo,y=score,fill=fone)) +
stat_pvalue_manual(tukeyfull,label="p.adj",y.position = c(4:11),size=2)
Interactions do not map onto plot as desired
Note: I know that it's advisable to do something like tukeyfull %>% add_xy_position(y="ftwo") before adding the test's p-values, but I get this error when I try to run that example: Error: Must group by variables found in '.data'. * Column 'fone:ftwo' is not found.
