专题三、方差分析
Wei Wei
2022-10-13
教学目标
- 掌握单因素方差分析、双因素方差分析、协方差分析方法、掌握重复测量方差分析方法
- 掌握混合设计的方差分析方法
方差分析
单因素方差分析
高考时填报学前教育志愿的顺序对于专业认同得分是否有影响？
因变量：专业认同得分
自变量：高考时填报学前教育志愿的顺序（三个水平）
第一步，进行方差齐性检验。
#导入数据
# 直接file-》import dataset，选择文件即可
library(haven)
# 导入数据默认为数据框格式（data.frame）
eee_data <- read_sav("https://gitee.com/vv_victorwei/r-language-data-analysis/raw/master/Untitled2.sav")
# head会呈现eee_data的前六行数据，以预览概况
head(eee_data)## # A tibble: 6 × 55
## V1 time_used IP id gender age birth…¹ roman…² hukou hair homet…³
## <chr> <dbl> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 175 594 223.66… 2225… 1 23 1 0 0 1 0
## 2 124 178 101.88… 2222… 1 23 1 0 0 1 0
## 3 138 207 211.16… 2225… 1 22 0 0 1 1 0
## 4 74 176 58.246… 2122… 1 23 1 0 1 1 0
## 5 149 232 120.25… 2222… 1 22 0 0 0 1 0
## 6 11 142 58.40.… 1601… 1 23 1 0 1 1 0
## # … with 44 more variables: family_income <dbl>, minzhu <dbl>, zhiyuan <dbl>,
## # height <dbl>, weight <dbl>, zhengzhi <dbl>, personality <dbl>, a1 <dbl>,
## # a2 <dbl>, a3 <dbl>, a4 <dbl>, a5 <dbl>, a6 <dbl>, a7 <dbl>, a8 <dbl>,
## # a9 <dbl>, a10 <dbl>, a11 <dbl>, a12 <dbl>, a13 <dbl>, b1 <dbl>, b2 <dbl>,
## # b3 <dbl>, b4 <dbl>, b5 <dbl>, b6 <dbl>, importance_ranking <chr>,
## # importance_ranking_d <dbl>, daode <dbl>, zhishi <dbl>, jineng <dbl>,
## # yishu <dbl>, food <chr>, chuan_food <dbl>, yue_food <dbl>, …# 计算专业认同所有题目的均值
eee_data$zyrt<-rowMeans(eee_data[19:31],na.rm = TRUE)
# 顺序变量转化为因子的函数ordered
eee_data$zhiyuan_f<-ordered(eee_data$zhiyuan,levels=c(1,2,3),labels = c('first choice','second choice','thrid choice'))
head(eee_data$zhiyuan_f)## [1] thrid choice thrid choice first choice second choice thrid choice
## [6] thrid choice
## Levels: first choice < second choice < thrid choice# 使用car库的拉文检验函数
library(car)## Loading required package: carDataleveneTest(zyrt ~ zhiyuan_f, data = eee_data)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 4.0165 0.01973 *
## 172
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1结果显著，说明方差不齐性
- 第二步，进行单因素方差分析。
# 单因素方差分析F检验结果,由于方差不齐性，var.equal = FALSE
oneway.test(zyrt~zhiyuan_f,data=eee_data,var.equal = FALSE)##
## One-way analysis of means (not assuming equal variances)
##
## data: zyrt and zhiyuan_f
## F = 0.49872, num df = 2.000, denom df = 83.444, p-value = 0.6091# 方差分析表
anova_zyrt<-aov(zyrt~zhiyuan_f,data=eee_data)
summary(anova_zyrt)## Df Sum Sq Mean Sq F value Pr(>F)
## zhiyuan_f 2 0.18 0.09065 0.487 0.615
## Residuals 172 32.02 0.18619# 计算效应量 eta square
library(lsr)
etaSquared(anova_zyrt)## eta.sq eta.sq.part
## zhiyuan_f 0.005629165 0.005629165# 绘制图表
boxplot(zyrt~zhiyuan_f,data=eee_data)
结果不显著，F = 0.49872, num df = 2.000, denom df = 83.444, p-value = 0.6091，说明三组的专业认同没有差异。
课堂例子
library(haven)
e1_3groups <- read_sav("https://gitee.com/vv_victorwei/r-language-data-analysis/raw/master/%E5%9D%87%E5%80%BC%E5%88%86%E6%9E%90/4.6%20E1%20three%20groups%20comparison.sav")
head(e1_3groups)## # A tibble: 6 × 2
## grade testscore
## <dbl+lbl> <dbl>
## 1 1 [小班] 70
## 2 1 [小班] 80
## 3 1 [小班] 75
## 4 1 [小班] 86
## 5 1 [小班] 77
## 6 1 [小班] 75e1_3groups$grade_f<-factor(e1_3groups$grade,levels = c(1,2,3),labels = c("class 1","class 2","class 3"))
library(car)
leveneTest(testscore ~ grade_f, data = e1_3groups)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 0.5671 0.5756
## 21# 单因素方差分析F检验结果,由于方差齐性，var.equal = TRUE
oneway.test(testscore ~ grade_f, data = e1_3groups,var.equal = TRUE)##
## One-way analysis of means
##
## data: testscore and grade_f
## F = 6.4131, num df = 2, denom df = 21, p-value = 0.006701# 方差分析表
anova_class<-aov(testscore ~ grade_f, data = e1_3groups)
summary(anova_class)## Df Sum Sq Mean Sq F value Pr(>F)
## grade_f 2 618.2 309.1 6.413 0.0067 **
## Residuals 21 1012.2 48.2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 计算效应量 eta square
library(lsr)
etaSquared(anova_class)## eta.sq eta.sq.part
## grade_f 0.3791782 0.3791782# 绘制图表
boxplot(testscore ~ grade_f, data = e1_3groups)
#事后检验
TukeyHSD(anova_class)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = testscore ~ grade_f, data = e1_3groups)
##
## $grade_f
## diff lwr upr p adj
## class 2-class 1 2.750 -5.9998993 11.4999 0.7118073
## class 3-class 1 11.875 3.1251007 20.6249 0.0069489
## class 3-class 2 9.125 0.3751007 17.8749 0.0399619#图示事后检验,虚线如果在置信区间之外，则结果显著
plot(TukeyHSD(anova_class))
多因素方差分析
分析地域（南方，北方）和报告志愿（第一，第二，第三）是否对专业认同有影响。 自变量：地域和志愿 因变量：专业认同
# 导入数据
# 直接file-》import dataset，选择文件即可
library(haven)
# 导入数据默认为数据框格式（data.frame）
eee_data <- read_sav("https://gitee.com/vv_victorwei/r-language-data-analysis/raw/master/Untitled2.sav")
# head会呈现eee_data的前六行数据，以预览概况
head(eee_data)## # A tibble: 6 × 55
## V1 time_used IP id gender age birth…¹ roman…² hukou hair homet…³
## <chr> <dbl> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 175 594 223.66… 2225… 1 23 1 0 0 1 0
## 2 124 178 101.88… 2222… 1 23 1 0 0 1 0
## 3 138 207 211.16… 2225… 1 22 0 0 1 1 0
## 4 74 176 58.246… 2122… 1 23 1 0 1 1 0
## 5 149 232 120.25… 2222… 1 22 0 0 0 1 0
## 6 11 142 58.40.… 1601… 1 23 1 0 1 1 0
## # … with 44 more variables: family_income <dbl>, minzhu <dbl>, zhiyuan <dbl>,
## # height <dbl>, weight <dbl>, zhengzhi <dbl>, personality <dbl>, a1 <dbl>,
## # a2 <dbl>, a3 <dbl>, a4 <dbl>, a5 <dbl>, a6 <dbl>, a7 <dbl>, a8 <dbl>,
## # a9 <dbl>, a10 <dbl>, a11 <dbl>, a12 <dbl>, a13 <dbl>, b1 <dbl>, b2 <dbl>,
## # b3 <dbl>, b4 <dbl>, b5 <dbl>, b6 <dbl>, importance_ranking <chr>,
## # importance_ranking_d <dbl>, daode <dbl>, zhishi <dbl>, jineng <dbl>,
## # yishu <dbl>, food <chr>, chuan_food <dbl>, yue_food <dbl>, …# 对hometown进行因子化，0南方，1北方
eee_data$hometown_f<-factor(eee_data$hometown,levels=c(0,1),labels = c('south','north'))
head(eee_data$hometown)## [1] 0 0 0 0 0 0head(eee_data$hometown_f)## [1] south south south south south south
## Levels: south north# 顺序变量转化为因子的函数ordered
eee_data$zhiyuan_f<-ordered(eee_data$zhiyuan,levels=c(1,2,3),labels = c('first choice','second choice','thrid choice'))
head(eee_data$zhiyuan)## [1] 3 3 1 2 3 3head(eee_data$zhiyuan_f)## [1] thrid choice thrid choice first choice second choice thrid choice
## [6] thrid choice
## Levels: first choice < second choice < thrid choice# 对每一行，计算专业认同所有题目均值
eee_data$zyrt<-rowMeans(eee_data[19:31],na.rm = TRUE)
# 第一步，方差齐性检验
# zyrt~hometown_f*zhiyuan_f
# ‘+’表示前后两个自变量/因素
# ‘:’表示墙后两个自变量/因素的交互作用
# ‘*’表示两个自变量/因素+二者的交互作用
# zyrt~hometown_f * zhiyuan_f表示zyrt是因变量，hometown_f和zhiyuan_f是自变量，*表示hometown_f+zhiyuan_f+hometown_f:zhiyuan_f
library(car)
leveneTest(zyrt~hometown_f*zhiyuan_f, data=eee_data)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 5 2.312 0.0461 *
## 169
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 第二步，方差分析
anova2_zyrt<-aov(zyrt~hometown_f*zhiyuan_f,data=eee_data)
summary(anova2_zyrt)## Df Sum Sq Mean Sq F value Pr(>F)
## hometown_f 1 0.05 0.04579 0.244 0.622
## zhiyuan_f 2 0.17 0.08701 0.464 0.630
## hometown_f:zhiyuan_f 2 0.30 0.14833 0.791 0.455
## Residuals 169 31.69 0.18751# 计算效应量eta squared
library(lsr)
etaSquared(anova2_zyrt)## eta.sq eta.sq.part
## hometown_f 0.001196169 0.001214188
## zhiyuan_f 0.005403622 0.005461696
## hometown_f:zhiyuan_f 0.009211167 0.009274468#绘图
boxplot(zyrt~hometown_f*zhiyuan_f,data=eee_data)
#绘制主效应和交互作用图
library(HH)## Loading required package: lattice## Loading required package: grid## Loading required package: latticeExtra## Loading required package: multcomp## Loading required package: mvtnorm## Loading required package: survival## Loading required package: TH.data## Loading required package: MASS##
## Attaching package: 'TH.data'## The following object is masked from 'package:MASS':
##
## geyser## Loading required package: gridExtra##
## Attaching package: 'HH'## The following objects are masked from 'package:car':
##
## logit, vif
interaction2wt(zyrt~hometown_f*zhiyuan_f,data=eee_data)
## 课堂例子
library(haven)
e2_2way <- read_sav("https://gitee.com/vv_victorwei/r-language-data-analysis/raw/master/%E5%9D%87%E5%80%BC%E5%88%86%E6%9E%90/two%20factors%20ANOVA(1).sav")
head(e2_2way)## # A tibble: 6 × 4
## id gender TV score
## <dbl> <dbl+lbl> <dbl+lbl> <dbl>
## 1 1 1 [girl] 1 [less than 15mins] 60
## 2 2 1 [girl] 1 [less than 15mins] 75
## 3 3 1 [girl] 1 [less than 15mins] 55
## 4 4 1 [girl] 1 [less than 15mins] 50
## 5 5 1 [girl] 1 [less than 15mins] 60
## 6 6 1 [girl] 1 [less than 15mins] 75e2_2way$gender_f<-factor(e2_2way$gender,levels = c(1,2),labels = c('girl','boy'))
e2_2way$TV_f<-factor(e2_2way$TV,levels = c(1,2,3),labels = c('<15 mins','15-30 mins','>30 mins'))
# 第一步，方差齐性检验
library(car)
leveneTest(score~TV_f*gender_f, data=e2_2way)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 5 0.6569 0.658
## 42# 第二步，方差分析
anova2<-aov(score~TV_f*gender_f, data=e2_2way)
summary(anova2)## Df Sum Sq Mean Sq F value Pr(>F)
## TV_f 2 1004 502.1 4.005 0.0256 *
## gender_f 1 13 13.0 0.104 0.7488
## TV_f:gender_f 2 879 439.6 3.506 0.0391 *
## Residuals 42 5266 125.4
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 计算效应量eta squared
library(lsr)
etaSquared(anova2)## eta.sq eta.sq.part
## TV_f 0.14020798 0.1601595
## gender_f 0.00181805 0.0024667
## TV_f:gender_f 0.12275471 0.1430751#简单效应分析
library(emmeans)##
## Attaching package: 'emmeans'## The following object is masked from 'package:HH':
##
## as.glhtSimple.Effects.by.TV<-emmeans(anova2, ~gender_f|TV_f)
pairs(Simple.Effects.by.TV,adjust='tukey')## TV_f = <15 mins:
## contrast estimate SE df t.ratio p.value
## girl - boy -5.62 5.6 42 -1.005 0.3208
##
## TV_f = 15-30 mins:
## contrast estimate SE df t.ratio p.value
## girl - boy -4.38 5.6 42 -0.781 0.4389
##
## TV_f = >30 mins:
## contrast estimate SE df t.ratio p.value
## girl - boy 13.12 5.6 42 2.344 0.0239Simple.Effects.by.gender<-emmeans(anova2, ~TV_f|gender_f)
pairs(Simple.Effects.by.gender,adjust='tukey')## gender_f = girl:
## contrast estimate SE df t.ratio p.value
## <15 mins - (15-30 mins) -1.25 5.6 42 -0.223 0.9729
## <15 mins - >30 mins 0.00 5.6 42 0.000 1.0000
## (15-30 mins) - >30 mins 1.25 5.6 42 0.223 0.9729
##
## gender_f = boy:
## contrast estimate SE df t.ratio p.value
## <15 mins - (15-30 mins) 0.00 5.6 42 0.000 1.0000
## <15 mins - >30 mins 18.75 5.6 42 3.349 0.0048
## (15-30 mins) - >30 mins 18.75 5.6 42 3.349 0.0048
##
## P value adjustment: tukey method for comparing a family of 3 estimates#绘图
boxplot(score~TV_f*gender_f, data=e2_2way)
#交互作用图
interaction.plot(x.factor = e2_2way$TV_f, trace.factor = e2_2way$gender_f,
response = e2_2way$score, fun = mean,
type = "b", legend = TRUE,
xlab = "TV times", ylab="Score",
pch=c(1,19), col = c("#00AFBB", "#E7B800"))
#绘制主效应和交互作用图
library(HH)
interaction2wt(score~TV_f*gender_f, data=e2_2way)
协方差分析
- 协方差分析是在控制了无关变量的基础上，考察自变量的效应。做法很简单，就是双因素方差分析而已
- 协变量可以是连续变量，也可以是分类变量
- 区别1：没有交互作用
- 区别2：协变量放在第一个
- 因此代码是weight~height_corrected + hometown_f, data=eee_data
- 控制了身高之后，南方和北方学生的体重是否有差异？
library(haven)
# 导入数据默认为数据框格式（data.frame）
eee_data <- read_sav("https://gitee.com/vv_victorwei/r-language-data-analysis/raw/master/Untitled2.sav")
# head会呈现eee_data的前六行数据，以预览概况
head(eee_data)## # A tibble: 6 × 55
## V1 time_used IP id gender age birth…¹ roman…² hukou hair homet…³
## <chr> <dbl> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 175 594 223.66… 2225… 1 23 1 0 0 1 0
## 2 124 178 101.88… 2222… 1 23 1 0 0 1 0
## 3 138 207 211.16… 2225… 1 22 0 0 1 1 0
## 4 74 176 58.246… 2122… 1 23 1 0 1 1 0
## 5 149 232 120.25… 2222… 1 22 0 0 0 1 0
## 6 11 142 58.40.… 1601… 1 23 1 0 1 1 0
## # … with 44 more variables: family_income <dbl>, minzhu <dbl>, zhiyuan <dbl>,
## # height <dbl>, weight <dbl>, zhengzhi <dbl>, personality <dbl>, a1 <dbl>,
## # a2 <dbl>, a3 <dbl>, a4 <dbl>, a5 <dbl>, a6 <dbl>, a7 <dbl>, a8 <dbl>,
## # a9 <dbl>, a10 <dbl>, a11 <dbl>, a12 <dbl>, a13 <dbl>, b1 <dbl>, b2 <dbl>,
## # b3 <dbl>, b4 <dbl>, b5 <dbl>, b6 <dbl>, importance_ranking <chr>,
## # importance_ranking_d <dbl>, daode <dbl>, zhishi <dbl>, jineng <dbl>,
## # yishu <dbl>, food <chr>, chuan_food <dbl>, yue_food <dbl>, …# 因子化hometown
eee_data$hometown_f = factor(eee_data$hometown,levels = c(0,1),labels = c('south','north'))
# 矫正原始数据的身高值
eee_data$height_corrected<-ifelse(eee_data$height>100,0.01*eee_data$height,eee_data$height)
head(eee_data$height_corrected)## [1] 1.62 1.67 1.64 1.68 1.68 1.66# 第一步，方差齐性检验
library(car)
leveneTest(weight~hometown_f, data=eee_data)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.5366 0.4648
## 173# 第二步，方差分析
anova2<-aov(weight~height_corrected + hometown_f, data=eee_data)
summary(anova2)## Df Sum Sq Mean Sq F value Pr(>F)
## height_corrected 1 4265 4265 74.69 3.7e-15 ***
## hometown_f 1 1 1 0.01 0.92
## Residuals 172 9821 57
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 计算效应量eta squared
library(lsr)
etaSquared(anova2)## eta.sq eta.sq.part
## height_corrected 3.026693e-01 3.027057e-01
## hometown_f 4.103655e-05 5.885474e-05重复测量方差分析
- 对于两组相关样本，使用配对样本t检验；对于三组相关样本，需要采用重复测量方差分析
library(haven)
# 导入数据默认为数据框格式（data.frame）
three_waves <- read_sav("https://gitee.com/vv_victorwei/r-language-data-analysis/raw/master/%E5%9D%87%E5%80%BC%E5%88%86%E6%9E%90/4.6%20E3%20R%20ANOVA.sav")
# head会呈现的前六行数据，以预览概况
head(three_waves)## # A tibble: 6 × 3
## score_t1 score_t2 score_t3
## <dbl> <dbl> <dbl>
## 1 64 75 76
## 2 62 70 71
## 3 70 72 68
## 4 60 73 74
## 5 89 92 93
## 6 76 82 83# 重复测量的数据变换
# R需要的模式是ID,TIMEPOINT, SCORE
# 可以把t1、t2和t3的数据堆砌一下，ID按1：20重复3遍堆砌，timepoint按1重复20遍、2重复20次，3重复20次堆砌
id<-rep(1:20,3)
time<-c(rep(1,20),rep(2,20),rep(3,20))
score<-c(three_waves$score_t1,three_waves$score_t2,three_waves$score_t3)
# 因子化
time<-factor(time)
id<-factor(id)
#合并三个变量为一个数据框
three_waves_new<-data.frame(id,time,score)
# 描述统计
library(psych)##
## Attaching package: 'psych'## The following object is masked from 'package:HH':
##
## logit## The following object is masked from 'package:car':
##
## logitdescribeBy(three_waves_new$score,group=three_waves_new$time)##
## Descriptive statistics by group
## group: 1
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 20 74.85 9.8 74.5 74.44 10.38 60 92 32 0.23 -1.13 2.19
## ------------------------------------------------------------
## group: 2
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 20 78.8 7.74 76.5 78.06 7.41 70 95 25 0.59 -0.96 1.73
## ------------------------------------------------------------
## group: 3
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 20 79.2 8.19 76.5 78.56 8.15 68 96 28 0.53 -0.98 1.83#图示发展变化
library(ggplot2)##
## Attaching package: 'ggplot2'## The following objects are masked from 'package:psych':
##
## %+%, alpha## The following object is masked from 'package:latticeExtra':
##
## layerggplot(three_waves_new, aes(x = time, y = score)) + geom_line(aes(group = id))
# 方差分析
library(ez)
three_waves.model <- ezANOVA(data = three_waves_new, dv = score, wid = id, within = time)
three_waves.model## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 time 2 38 16.26595 7.88065e-06 * 0.05175108
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 2 time 0.2696767 7.543823e-06 *
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
## 2 time 0.5779267 0.000349976 * 0.5917085 0.0003088779 *library(rstatix)##
## Attaching package: 'rstatix'## The following object is masked from 'package:MASS':
##
## select## The following object is masked from 'package:stats':
##
## filterpairwise_t_test(data = three_waves_new,score ~ time, paired = TRUE,p.adjust.method = 'bonferroni')## # A tibble: 3 × 10
## .y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 score 1 2 20 20 -4.12 19 0.00058 0.002 **
## 2 score 1 3 20 20 -4.14 19 0.00055 0.002 **
## 3 score 2 3 20 20 -1.19 19 0.248 0.744 ns混合设计的方差分析
- 研究设计同时包含被试间设计和被试内设计
# 实验组对照组接受前后测，组间（实验vs.对照），组内（前vs.后测）
library(haven)
# 导入数据默认为数据框格式（data.frame）
rct_data <- read_sav("https://gitee.com/vv_victorwei/r-language-data-analysis/raw/master/%E5%9D%87%E5%80%BC%E5%88%86%E6%9E%90/4.6%20R4%20RCT%20data.sav")
# head会呈现的前六行数据，以预览概况
head(rct_data)## # A tibble: 6 × 5
## age group pretest posttest grade
## <dbl> <dbl+lbl> <dbl> <dbl> <dbl>
## 1 50 0 [control] 35 35 1
## 2 51 0 [control] 34 35 1
## 3 51 0 [control] 26 30 1
## 4 53 0 [control] 31 32 1
## 5 55 0 [control] 36 34 1
## 6 56 0 [control] 36 36 1#转换数据类型，score只要一列
score<-c(rct_data$pretest,rct_data$posttest)
id<-factor(rep(1:32,2))
time<-factor(c(rep(0,32),rep(1,32)))
group<-factor(rep(rct_data$group,2))
grade<-factor(rep(rct_data$grade,2))
#变量合并为一个数据框
rct_new<-data.frame(score,id,group,time,grade)
library(ez)
rct_model<-ezANOVA(data = rct_new, dv = score, wid = id,
between = group, within = .(time), type = 3, detailed = F)
rct_model## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 group 1 30 4.27111 4.749097e-02 * 0.11387829
## 3 time 1 30 68.75250 2.959674e-09 * 0.18237895
## 4 group:time 1 30 18.45309 1.684530e-04 * 0.05648734#事后检验/简单效应分析
#编码组别：1 对照组-前测，2对照组-后测，3实验组-前测，4实验组-后测
rct_new$groups[rct_new$group==0 & rct_new$time==0]<-1
rct_new$groups[rct_new$group==0 & rct_new$time==1]<-2
rct_new$groups[rct_new$group==1 & rct_new$time==0]<-3
rct_new$groups[rct_new$group==1 & rct_new$time==1]<-4
library(rstatix)
pairwise.t.test(rct_new$score, rct_new$groups, paired = TRUE, p.adjust.method = "bonferroni")##
## Pairwise comparisons using paired t tests
##
## data: rct_new$score and rct_new$groups
##
## 1 2 3
## 2 0.0159 - -
## 3 1.0000 1.0000 -
## 4 0.0044 0.0282 1e-05
##
## P value adjustment method: bonferroni