【R语言】方差分析

#载入示例数据
data(iris)

#载入包
library(tidyverse) 
library(lme4)  

#生成个体id
iris$ID <- rep(1:30,each=5)  #假设一朵花测五次
iris$ID <- as.factor(iris$ID)

#-------------【ANOVA】------------
iris_aov1<-aov(Sepal.Length~Species,data=iris)
summary(iris_aov1)

iris_aov2<-aov(Sepal.Length~Species+Error(ID),data=iris)
summary(iris_aov2)

iris_aov3<-lm(Sepal.Length~Species,data=iris)
summary(iris_aov3) #R2又称为方程的确定性系数（coefficient of determination）
anova(iris_aov3) 
anova(iris_aov3) ;summary(iris_aov1) #结果同 iris_aov1

#考虑ID为随机效应

library(lme4)
#固定
iris_lmer1 <- lmer(Sepal.Length~Species+(1|ID),data = iris)
#每一个ID在不同Species的截距不同
iris_lmer2 <- lmer(Sepal.Length~Species+(1+(Species|ID)),data = iris)
summary(iris_lmer1)
summary(iris_lmer2)

library(sjPlot)
coef(iris_lmer1)
plot(ranef(iris_lmer1))  # 画随机效应
plot(iris_lmer1) #residuals残差
plot_model(iris_lmer1, type = "re", show.values = TRUE)
#查看每个ID分别计算的截距和斜率
coefficients(iris_lmer1)

# plot the residuals, it is good that all data points are around 0 and symetric, 
#whith no obvious pattern
#The random effect part tells you how much variance you find 
#   among levels of our grouping factors(s), plus the residual variance.
#0.261/(0.261+0.004)=??? (The differences between speakers 
#        explain ???% of the variance that’s “left over”
#         after the variance explained by the fixed effects.)

#T value: Estimate/SE   T>1.96 indicates that p<0.05




# 比较模型与去除固定效应模型的差异，若2个模型差异显著，则认为固定效应存在。
# 采用 maximun likelihood去比较
iris_lmer1_2 <- lmerTest::lmer(REML = F,Sepal.Length~Species+(1|ID),iris)
iris_lmer2_2 <- lmer(REML = F,Sepal.Length~Species+(1+(Species|ID)),iris)
anova(iris_lmer1_2,iris_lmer2_2)


#post-hoc test in lmeTest 事后两两比较
iris_post <- emmeans::emmeans(iris_lmer1,pairwise~Species,
                              adjust="none")
iris_post


#【探索性作图】
#作图查看分布情况
plot_iris <- ggplot(data = iris, 
                    aes(x = Species, y = Sepal.Length,fill=Species))    %>%
  +geom_boxplot(position=position_dodge(0.8),size=0.1)  %>%
  +labs(x =element_blank(), y = "Sepal.Length")+theme_classic()

plot_iris

#
# 散点图
plot_iris2 <- ggplot(data = iris, 
                     aes(x = Species, y = Sepal.Length,color=Species))    %>%
  +geom_dotplot(aes(color = Species,fill=Species),
                stackdir='center', #散点中心对称
                binaxis = "y", #binaxis="y"是指沿着y轴进行分箱
                dotsize=0.8)  %>%
  +labs(x =element_blank(),y = "Sepal.Length")+theme_classic()
plot_iris2

#箱线图
plot_iris+ geom_dotplot(stackdir='center', binaxis = "y",dotsize=0.8)

#重新排序变量的分类
#reorder the levels as the way you wish, otherwise it is alphabet order
iris$Species=factor(iris$Species,levels=c("versicolor","virginica","setosa"))


# 绘制所有数据，查看是否存在任何极值
plot(iris$Sepal.Length)
# 从最小值到最大值显示数据
plot(sort(iris$Sepal.Length),ylab="Sepal.Length)") 

# 检查数据是否符合正态分布
hist(iris$Sepal.Length) # Plot the histrogram figure
qqnorm(iris[iris$Species=="versicolor", ]$Sepal.Length)
qqline(iris[iris$Species=="versicolor", ]$Sepal.Length,col="red")



##作图 展示主效应和交互效应
plot_Species <- ggplot(data = iris, aes(x = Species, y = Sepal.Length)) %>%
     +geom_boxplot(position=position_dodge(0.8),size=0.1)  %>%
     +labs(x =element_blank(), y = "Sepal.Length")+theme_classic()
plot1


plot_ID <- ggplot(data = iris, aes(x = ID, y = Sepal.Length))   %>%
          +geom_boxplot(aes(fill = NULL),
                        position=position_dodge(0.8),size=0.1)  %>%
          + scale_color_brewer(palette="Dark2")   %>%
          +labs(x =element_blank(),y = "Sepal.Length")+theme_classic()
plot_ID

plot_interaction<-ggplot(data = iris, aes(x = Species, y = Sepal.Length)) %>%
     +geom_boxplot(aes(fill = ID),position=position_dodge(0.8),size=0.1) %>%
     + scale_color_brewer(palette="Dark2")+labs(x =element_blank(),
                                        y = "Sepal.Length")+theme_classic()
plot_interaction


plot_interaction2 <- ggplot(iris,aes(x=Species, y = Sepal.Length,
                                     group=ID,color=ID))  %>%
                     +geom_line()+geom_point()
plot_interaction2