image.png image.png

library(tidyverse)
library(ggplot2)
library(ggsci)
library(ggsignif)
library(ggpubr)
library(car)
library(devtools)
library(userfriendlyscience)


data("ToothGrowth")
 ToothGrowth$dosf <- factor(x = ToothGrowth$dose, ordered = T)
  

1.单因素方差分析

1.1 正态性检验

 with(ToothGrowth, tapply(len,dosf,shapiro.test))

1.2 方差齐性检验

  leveneTest(len~dosf, ToothGrowth)

1.3 单因素ANOVA

  AOV1 <- aov(len~dosf, ToothGrowth)
  summary(AOV1)

1.4 诊断模型（残差检验）

  res1 <- residuals(AOV1)
  shapiro.test(res1)  #指标：差异不显著
  leveneTest(res1~dosf,ToothGrowth) #指标：方差齐
  ggqqplot(res1)
  

1.5两两比较

  Tu <- TukeyHSD(AOV1)
  P <- Tu$dosf[,4] 
  format(P, scientific = T)
  
  

1.6 作图（不要用ggboxplot）

### ggboxplot(
#    ToothGrowth, x="dosf",y="len", fill="dosf",
#    palette = "npg" 
#  ) +
#  stat_compare_means(method = "anova")  ##不可用：P值未经过矫正，两两比较的结果
  

 ggplot(ToothGrowth,aes(x=dosf,y=len,fill=dosf)) +
  stat_boxplot(geom = "errorbar") +
  geom_boxplot()+
  geom_point() +
  geom_signif(comparisons = list(c("0.5","1"),c("0.5","2"),c("1","2")),
              y_position = c(36,40,38), #P值坐标
              annotations = format(unname(P),scientific = T, digits = 3))+
  #矫正后的P值；需与list的组对应
    labs(x=NULL,y="Tooth Length (cm)", title = "Tooth Growth", fill="Dose (mg/d)")  +
  scale_fill_npg() +
  theme_classic2()
              

2. Welch ANOVA(正态但方差不齐)

  oneway.test(len~dosf, ToothGrowth,var.equal = F)
  #var.equal=T 则为aov()
  
  ##两两比较
  GH <- with(ToothGrowth,posthocTGH(len,dosf, method="games-howel",digits=3))
  GH$output

3.正态不齐

  kruskal.test(len~dosf,ToothGrowth)
  ##两两比较
  with(ToothGrowth, pairwise.wilcox.test(len,dosf,p.adjust.method = "BH"))
  

4.双因素方差分析

  AOV2 <- aov(len~dosf*supp,ToothGrowth) # *表示二者有交互作用；等价于aov(len~dosf+supp+dosf:supp,ToothGrowth)
  summary(AOV2) #先看交互作用的显著性
  
  ##两两比较
 TK <-  TukeyHSD(AOV2,'dosf:supp')
 TKint <- as.data.frame(TK$`dosf:supp`) 
 TKint$pair <- rownames(TKint)
 TKint$sig <- cut(TKint$`p adj`,c(0,0.05,1),labels = c("P<0.05","NS"),right = F) #NS:not signaficant, right闭区间
 
 
  ##作图
  ggplot(ToothGrowth,aes(x=dosf,y=len,fill=supp)) +
  stat_boxplot(geom = "errorbar") +
  geom_boxplot()  #组内箱线图
  
  plot(TK) #森林图
  
  ggplot(TKint, aes(x=`p adj`, y=pair,color=sig)) +  ##p adj中间有空格，所以用``
  geom_errorbar(aes(xmin=lwr, xmax=upr),high=0.3) +
    geom_vline(xintercept = 0,lty=2) +
    labs(x=NULL,y=NULL, title = "95% family-wise confidence level",color="significant")  +
  scale_color_npg() +
    theme_classic2()