使用本地数据 鸢尾花（yuān wěi huā）初探线性混合模型Linear Mixde Model
施工ing 交互作用的
#载入示例数据
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_lmer) #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
|
