dta <- read.csv("rt_dummy_data.csv")
head(dta,20)
## PID RT modality stim
## 1 301 1024 Audio-only gown
## 2 301 838 Audio-only might
## 3 301 1060 Audio-only fern
## 4 301 882 Audio-only vane
## 5 301 971 Audio-only pup
## 6 301 1064 Audio-only rise
## 7 301 1283 Audio-only hill
## 8 301 982 Audio-only caught
## 9 301 1144 Audio-only leg
## 10 301 743 Audio-only hurt
## 11 301 639 Audio-only shack
## 12 301 853 Audio-only tell
## 13 301 649 Audio-only worth
## 14 301 757 Audio-only shook
## 15 301 958 Audio-only douse
## 16 301 942 Audio-only jail
## 17 301 748 Audio-only bog
## 18 301 994 Audio-only thatch
## 19 301 810 Audio-only bit
## 20 301 1241 Audio-only pop
library(pacman)
## Warning: 套件 'pacman' 是用 R 版本 4.2.3 來建造的
p_load(lme4,tidyverse,afex)
dta_anova <- aggregate(dta$RT,list(dta$PID,dta$modality),mean)
names(dta_anova) <- c('PID','modality','RTmean')
dta_anova$PID <- as.factor(dta_anova$PID)
head(dta_anova,10)
## PID modality RTmean
## 1 301 Audio-only 1027.4106
## 2 302 Audio-only 1047.3455
## 3 303 Audio-only 882.6735
## 4 304 Audio-only 1238.2754
## 5 306 Audio-only 1043.3231
## 6 307 Audio-only 1115.6829
## 7 308 Audio-only 1252.6053
## 8 309 Audio-only 795.7773
## 9 310 Audio-only 1175.5823
## 10 311 Audio-only 1015.0566
rst_anova <- aov(formula = RTmean ~ modality + Error(PID/modality),data=dta_anova)
summary(rst_anova)
##
## Error: PID
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 52 2928806 56323
##
## Error: PID:modality
## Df Sum Sq Mean Sq F value Pr(>F)
## modality 1 183264 183264 43.72 2.05e-08 ***
## Residuals 52 217966 4192
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
MNS <- tapply(dta_anova$RTmean,dta_anova$modality,mean)
MNS
## Audio-only Audiovisual
## 1044.073 1127.233
-兩組差 83 ms
p_load(effectsize)
eta_squared(rst_anova, partial = FALSE)
## # Effect Size for ANOVA (Type I)
##
## Group | Parameter | Eta2 | 95% CI
## ----------------------------------------------
## PID:modality | modality | 0.06 | [0.00, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
dta$modality <- ifelse(dta$modality == "Audio-only", 0, 1)
rst <- lmer(RT ~ 1 + modality + (1 + modality|PID) + (1 + modality|stim),
data = dta,
control = lmerControl(optimizer = "bobyqa"))
summary(rst)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: RT ~ 1 + modality + (1 + modality | PID) + (1 + modality | stim)
## Data: dta
## Control: lmerControl(optimizer = "bobyqa")
##
## REML criterion at convergence: 302385.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3646 -0.6964 -0.0140 0.5886 5.0003
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## stim (Intercept) 303.9 17.43
## modality 216.7 14.72 0.16
## PID (Intercept) 28552.7 168.98
## modality 7709.8 87.81 -0.17
## Residual 65258.8 255.46
## Number of obs: 21679, groups: stim, 543; PID, 53
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1044.14 23.36 52.14 44.704 < 2e-16 ***
## modality 83.18 12.58 52.09 6.615 2.02e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## modality -0.178