前言
- Two-Way ANOVA由于有两个自变量，在考虑两个变量单独对应变量影响的同时，还需要考虑交互作用
- 当ANOVA结果显示交互作用差异显著时，需要做每个变量在另一个变量各个水平下的显著性分析，而不能做变量的主效应分析（即使ANOVA结果显示单个变量存在显著性差异）
- 当交互作用差异不显著，而变量的差异显著时，才可以对相应的变量做主效应分析
R代码
######################################
## Two-Way Anova ##
####################################
library(psych)
library(car)
library(emmeans)
#Environment Clean
rm(list = ls())
#Import data
Input <- ("
Diet Country Weight_change
A USA 0.120
A USA 0.125
A USA 0.112
A UK 0.052
A UK 0.055
A UK 0.044
A NZ 0.080
A NZ 0.090
A NZ 0.075
B USA 0.096
B USA 0.100
B USA 0.089
B UK 0.025
B UK 0.029
B UK 0.019
B NZ 0.055
B NZ 0.065
B NZ 0.050
C USA 0.149
C USA 0.150
C USA 0.142
C UK 0.077
C UK 0.080
C UK 0.066
C NZ 0.055
C NZ 0.065
C NZ 0.050
C NZ 0.054
")
Data <- read.table(textConnection(Input),header=TRUE)
#Remove unnecessary objects
rm(Input)
#Check the data frame
headTail(Data) #from psych package
str(Data)
summary(Data)
tapply(Data[[3]],Data$Country ,mean) #看均值
#############################################
#运行结果
> headTail(Data) #from psych package
Diet Country Weight_change
1 A USA 0.12
2 A USA 0.12
3 A USA 0.11
4 A UK 0.05
... <NA> <NA> ...
25 C NZ 0.06
26 C NZ 0.06
27 C NZ 0.05
28 C NZ 0.05
> str(Data)
'data.frame': 28 obs. of 3 variables:
$ Diet : chr "A" "A" "A" "A" ...
$ Country : chr "USA" "USA" "USA" "UK" ...
$ Weight_change: num 0.12 0.125 0.112 0.052 0.055 0.044 0.08 0.09 0.075 0.096 ...
> summary(Data)
Diet Country Weight_change
Length:28 Length:28 Min. :0.01900
Class :character Class :character 1st Qu.:0.05350
Mode :character Mode :character Median :0.07050
Mean :0.07746
3rd Qu.:0.09700
Max. :0.15000
> tapply(Data[[3]],Data$Country ,mean)
NZ UK USA
0.06390000 0.04966667 0.12033333
#################################################
#判断是否存在交互，运行结果见图1
interaction.plot(x.factor = Data$Country,
trace.factor = Data$Diet,
response = Data$Weight_change,
fun = mean,
type="b",
col=c("black","red","green"), # Colors for levels of trace var.
pch=c(19, 17, 15), # Symbols for levels of trace var.
fixed=TRUE, # Order by factor order in data
leg.bty = "o")
#Specify the linear model and conduct an analysis of variance
model <- lm(Weight_change ~ Country + Diet + Country:Diet,
data = Data)
Anova(model, type = "II") #from car package
#################################################
#运行结果
Anova Table (Type II tests)
Response: Weight_change
Sum Sq Df F value Pr(>F)
Country 0.0256761 2 318.715 2.426e-15 ***
Diet 0.0051534 2 63.969 3.634e-09 ***
Country:Diet 0.0040162 4 24.926 2.477e-07 ***
Residuals 0.0007653 19
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##############################################
#Post-hoc testing with emmeans
#Main effects
#此例中交互的p值小于0.05，不应该做main effects，此处只是展示代码
em.resul1<-emmeans(model,
~ Diet, contr="tukey"); em.resul1
plot(em.resul1, comparisons = TRUE)
em.resul2<-emmeans(model,
~ Country, contr="tukey"); em.resul2
plot(em.resul2, comparisons = TRUE)
#############################################
#运行结果，可视化结果见图2
> #Post-hoc testing with emmeans
> #Main effects
> em.resul1<-emmeans(model,
+ ~ Diet, contr="tukey"); em.resul1
NOTE: Results may be misleading due to involvement in interactions
$emmeans
Diet emmean SE df lower.CL upper.CL
A 0.0837 0.00212 19 0.0792 0.0881
B 0.0587 0.00212 19 0.0542 0.0631
C 0.0924 0.00203 19 0.0882 0.0967
Results are averaged over the levels of: Country
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
A - B 0.02500 0.00299 19 8.356 <.0001
A - C -0.00878 0.00293 19 -2.997 0.0194
B - C -0.03378 0.00293 19 -11.533 <.0001
Results are averaged over the levels of: Country
P value adjustment: tukey method for comparing a family of 3 estimates
############################################
#interaction
em.resul3<-emmeans(model,
~ Country|Diet, contr="tukey"); em.resul3
plot(em.resul3, comparisons = TRUE)
em.resul4<-emmeans(model,
~ Diet|Country, contr="tukey"); em.resul4
plot(em.resul4, comparisons = TRUE)
##########################################
#运行结果，可视化结果见图3
> plot(em.resul2, comparisons = TRUE)
> #interaction
> em.resul3<-emmeans(model,
+ ~ Country|Diet, contr="tukey"); em.resul3
$emmeans
Diet = A:
Country emmean SE df lower.CL upper.CL
NZ 0.0817 0.00366 19 0.0740 0.0893
UK 0.0503 0.00366 19 0.0427 0.0580
USA 0.1190 0.00366 19 0.1113 0.1267
Diet = B:
Country emmean SE df lower.CL upper.CL
NZ 0.0567 0.00366 19 0.0490 0.0643
UK 0.0243 0.00366 19 0.0167 0.0320
USA 0.0950 0.00366 19 0.0873 0.1027
Diet = C:
Country emmean SE df lower.CL upper.CL
NZ 0.0560 0.00317 19 0.0494 0.0626
UK 0.0743 0.00366 19 0.0667 0.0820
USA 0.1470 0.00366 19 0.1393 0.1547
Confidence level used: 0.95
$contrasts
Diet = A:
contrast estimate SE df t.ratio p.value
NZ - UK 0.0313 0.00518 19 6.046 <.0001
NZ - USA -0.0373 0.00518 19 -7.204 <.0001
UK - USA -0.0687 0.00518 19 -13.251 <.0001
Diet = B:
contrast estimate SE df t.ratio p.value
NZ - UK 0.0323 0.00518 19 6.239 <.0001
NZ - USA -0.0383 0.00518 19 -7.397 <.0001
UK - USA -0.0707 0.00518 19 -13.637 <.0001
Diet = C:
contrast estimate SE df t.ratio p.value
NZ - UK -0.0183 0.00485 19 -3.782 0.0034
NZ - USA -0.0910 0.00485 19 -18.773 <.0001
UK - USA -0.0727 0.00518 19 -14.023 <.0001
P value adjustment: tukey method for comparing a family of 3 estimates
##########################################################
#p值可视化，结果见图4
em.resul5<-emmeans(model,
~ Diet:Country); em.resul5
pwpp(em.resul5,type = "response",by='Country')
图1：曲线有交叉说明存在交互，无交叉说明不存在交互。此例中，Country和Diet之间有交互。
图2：横坐标为均值，纵坐标为分组，蓝色条块代表各组的置信区间，红色箭头用于判断组间的显著性差异，当不同组的箭头重叠时，则不存在显著性差异。
图3
图4：横坐标为p值，纵坐标为分组。例如，country为NZ时，dietA和B、C间均存在显著性差异，且p值<0.001,而B和C无显著性差异，且p值接近1
